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Summary 


We  describe  a  theoretical  model  for  the  digital  simulation  of  flame¬ 
spreading  and  pressure  wave  propagation  in  axisymmetric  two  dimensional  pro¬ 
pelling  charges.  Because  longitudinal  pressure  waves  are  known  to  be  strong¬ 
ly  influenced  by  the  distribution  of  ullage  around  the  charge,  we  describe 
an  approach  which  recognizes  explicitly  all  internal  boundaries  defined  by 
jumps  in  porosity.  Likewise,  the  formulation  of  the  analysis  of  the  internal 
boundaries  enables  the  explicit  recognition  of  the  influence  on  flamespread¬ 
ing  of  bag  material,  liners  and  additives. 

An  important  milestone  in  respect  to  the  numerical  implementation  of 
this  model  is  the  development  of  a  technique  of  numerical  solution  for  a 
geometrically  complex  region  of  two  dimensional  two  phase  flow.  We  present 
a  solution  technique  based  on  an  equipotential  map  of  the  physical  domain 
onto  a  square  and  an  update  of  the  state  variables  by  an  explicit  two  step 
marching  technique.  A  characteristic  formulation  of  the  equations  is  used 
at  the  boundaries. 

Stable  solutions  are  presented  for  three  problems  based  on  nominal  data. 
All  three  problems  recognize  the  hemispherical  shape  of  the  breech  closure 
plug  and  the  taper  of  the  tube.  The  third  problem  considers  the  projectile 
to  have  a  boat tail  which  intrudes  into  the  propelling  charge. 
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1.0  INTRODUCTION 


This  report  presents  details  of  the  formulation  and  partial  development 
of  a  model  of  the  macroscopic,  two  dimensional,  two  phase,  unsteady,  react¬ 
ing  flow  occurring  in  a  gun. 

The  objectives  of  our  enquiry  and  the  scope  of  the  present  study  are 
described  in  section  1.1.  In  section  1.2  we  provide  some  background  infor¬ 
mation  which  coordinates  the  present  study  with  previous  work.  Finally,  in 
section  1.3  of  this  introduction,  we  summarize  our  approach  and  present 
findings . 

1 . 1  Objectives  and  Scope  of  Present  Study 

We  describe  first  of  all,  the  goals  of  our  research  as  a  whole.  Al¬ 
though  the  present  document  constitutes  a  final  contract  report,  it  is,  in 
a  sense,  an  interim  report-  For  this  reason  we  shall  try  to  clarify,  as 
much  as  possible,  the  extent  to  which  the  scope  of  the  present  study  permits 
us  to  reach  our  overall  goals  and  to  identify  those  aspects  of  the  model 
which  are  deferred  for  future  study. 

In  the  general  area  of  interior  ballistic  phenomena  in  medium  caliber 
guns  we  may  identify  three  topics  of  importance  to  the  charge  designer. 

These  are: 

(a)  The  central  problem  of  classical  interior  ballistics  ;  namely,  to  ob¬ 
tain  a  desired  muzzle  velocity  within  a  specified  pressure  limitation.  Of 
course,  the  designer  must  also  take  into  account  problems  of  systems  com¬ 
patibility  and  other  constraints,  but  we  are  concerned  here  only  with  the 
problems  which  have  some  hydrodynamic  content. 

(b)  The  minimization  of  tube  erosion  and  of  muzzle  blast. 

(c)  The  elimination  of  ignition  anomalies  which  appear,  in  their  mildest 
form,  as  longitudinal  pressure  waves  and,  in  their  worst  form,  as  cata¬ 
strophic  overpressures  of  the  gun. 

It  is  this  last  problem  which  is  addressed  by  our  enquiry  as  a  whole. 

Of  particular  interest  are  the  important  roles  played  by  the  detailed  dis¬ 
tribution  of  ullage,  the  igniter  discharge  characteristics  and  the  flow  in¬ 
hibition  due  to  the  presence  of  bag  materials  together  with  the  liners  and 
additives  used  to  control  erosion  and  wear. 

The  scope  of  the  present  study  is  as  follows.  We  provide  the  details 
of  a  model  of  the  macroscopic,  two  dimensional,  two  phase  unsteady  reacting 
flow  in  a  gun.  The  discussion  includes  details  of  the  governing  equations 
and  the  method  of  solution.  Numerical  solutions  are  generated  for  the  prob¬ 
lem  of  two  dimensional  convective  flamespreading  in  a  packed  bed  of  granular 
propellant  in  an  ullage  free,  but  geometrically  complex,  container  with 
stationary  impermeable  walls. 

In  its  present  state  of  numerical  implementation,  the  model  is  there¬ 
fore  useful  for  the  investigation  of  flamespreading  in  ullage-free  case  am¬ 
munition  in  which  the  projectile  intrusion  due  to  the  boattail  or  afterbody 
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may  be  significant.  However,  the  application  of  the  model  to  bag  charges, 
for  which  ullage  is  present,  and  to  the  sequence  of  events  following  the 
start  of  the  motion  of  the  projectile,  is  the  subject  of  future  work. 

1 . 2  Background 

The  ideal  propelling  charge  is  thought  to  be  one  in  which  ignition  oc¬ 
curs  simultaneously  at  all  points  so  that  all  exposed  surfaces  are  burning 
uniformly  from  the  initial  instant.  A  limiting  form  of  this  ideal  involves 
the  instantaneous  combustion  of  the  entire  charge.  From  a  theoretical  point 
of  view,  the  hydrodynamical  model  of  such  an  interior  ballistic  cycle  in¬ 
volves  the  study  of  a  single  phase  substance  expanding  from  a  quiescent  but 
energetic  initial  state.  Such  was  the  problem  considered  by  Lagrange  and 
later,  in  greater  detail,  by  Pidduck  and  Kent^ . 

From  a  practical  point  of  view,  however,  a  charge  which  burns  instan¬ 
taneously  is  one  which  subjects  the  tube  to  a  certain  level  of  stress  while 
propelling  the  projectile  at  much  lower  average  value.  Accordingly,  charges 
are  designed  to  release  their  energy  gradually,  producing  gas  at  a  rate  which 
compensates  for  the  motion  of  the  projectile.  In  this  fashion  system  *’cost” 
which  is  related  to  the  maximum  stress  is  better  aligned  with  system  ”effec- 
tiveness”  which  is  related  to  the  average  stress.  An  obvious  theoretical  con¬ 
sequence  of  this  procedure  is  that  a  properly  posed  model  of  interior  ballis¬ 
tic  phenomena  must  now  consider  the  coupled  motion  of  both  the  gas  and  solid 
phases.  By  treating  the  motion  as  that  of  a  well  stirred  mixture  subjected 
to  an  independently  determined  pressure  gradient  to  represent  the  response  to 
projectile  motion,  one  may  establish  a  lumped  parameter  model  of  interior 
ballistic  phenomena^. 

In  fact,  the  mixture  is  far  from  well  stirred  except  possibly  quite  late 
in  the  interior  ballistic  cycle  as  burnout  is  approached.  Since  many  modern 
charges  produce  burning  almost  throughout  the  entire  propulsion  event,  the 
assumption  of  well  stirredness  is  rather  limiting. 

In  order  to  understand  the  lack  of  equilibrium  between  the  phases  we  re¬ 
fer  to  figures  1.2.1  and  1.2.2  which  respectively  illustrate  a  typical  cased 
ammunition  charge  and  a  typical  bag  charge.  It  is  apparent  that  inhomogeneity 
exists  in  the  combustion  chamber  from  the  initial  instant  due  to  the  presence 
of  spaces  around  the  charge.  These  spaces  are  referred  to  as  ullage.  Then, 
when  the  igniter  begins  to  function,  further  inhomogeneities  arise  as  the 
charge  is  not  ignited  uniformly.  Rather,  a  spreading  process  occurs  which 
involves  a  front  referred  to  as  a  convective  flame.  Hot  pressurized  gas  pen¬ 
etrates  the  charge,  an  aggregate  of  regularly  formed  grains.  The  penetration 
is  accompanied  by  intense  convective  heat  transfer  which  ignites  the  charge, 
and  also  by  momentum  transfer  which  accelerates  the  charge.  Because  of  the 
impermeability  of  the  charge,  the  pressures  at  two  stations  in  the  chamber 
may  differ  by  tens  of  MPa  during  flamespreading. 

It  should  be  noted  that  the  ignition  system  depends  inherently  on  a  prin¬ 
ciple  of  inhomogeneous  venting  since  the  heat  transfer  rates  associated  with 
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linear  conduction  and  radiation  are  several  orders  of  magnitude  less  than 
those  associated  with  convection.  Yet,  the  pressure  gradients  associated 
with  inhomogeneous  venting  can  produce  pronounced  wave  phenomena  which  are 
superimposed  upon  the  ideal  pressurization  history.  Significantly,  the 
existence  of  longitudinal  pressure  waves  has  been  correlated  in  many  cases 
with  weapon  malfunction  due  to  overpressurization^ .  Therefore,  the  ignition 
train  is  frequently  designed  so  that,  as  shown  in  figures  1.2.1  and  1.2.2, 
venting  occurs  as  close  to  simultaneously  as  possible  over  the  length  of  the 
center  line.  The  inhomogeneity  is  therefore  concentrated  in  the  radial  di¬ 
rection  and  the  tendency  to  produce  longitudinal  waves  is  minimized. 


The  importance  of  igniter  location  and  of  the  distribution  of  ullage  has 
been  known  since  the  time  of  Kent^  and  of  Hedden  &  Nance^.  Kent  showed,  in 
the  context  of  bag  charges,  that  ullage  around  the  charge  was  most  influen¬ 
tial  in  reducing  the  overall  longitudinal  pressure  gradient  during  flame¬ 
spreading.  This  is  easy  to  understand  since  the  ullage  around  the  charge 
provides  a  low  impedance  path  for  the  flow  of  gas  from  the  high  pressure  sec¬ 
tions  to  the  low  pressure  sections.  Hedden  and  Nance  demonstrated  the  strong 
influence  exerted  by  longitudinal  ullage  on  wave  amplitude.  More  recently, 
several  papers  have  appeared  in  which  the  dependence  of  the  pressure  wave 
amplitude  on  primer /propellant  interface^ ’ ^ ,  charge  permeability®’®  and 
charge  configuration  or  ullage  distribution^ > 10  has  been  explored. 
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More  or  less  concurrently  with  these  experimental  studies,  the  modeling 
of  interior  ballistic  phenomena  according  to  an  unsteady,  two-phase,  contin¬ 
uum  approach  has  been  developed  to  the  point  of  being  able,  in  certain  cases, 
to  predict  the  structure  and  amplitude  of  the  longitudinal  pressure  waves. 
While  relatively  limited  flamespreading  models  were  advanced  by  several 
authors^^“^^ ,  the  models  of  Gough^^”^^  and  of  Fisher^^"^^  have  reflected  in¬ 
creasing  attention  to  the  detailed  configuration  of  the  charge. 

14-1  8 

The  NOVA  code  developed  by  Gough  effects  a  quasi-one-dimensional 

representation  of  the  charge  with  precise  attention  to  the  axial  distribution 
of  ullage.  This  code  has  shown  good  correlation  with  observatioyg  of  case 
ammunition^  5  as  well  as  with  certain  simpler  configurations'^  ,  The 

CALSPAN  code  developed  by  Fisher  goes  beyond  the  NOVA  code  in  terms  of 

phenomenological  completeness.  A  two  dimensional  capacity  exists,  although 
the  code  is  customarily  run  as  a  quasi-two-dimensional  model  in  which  one 
dimensional  flow  is  considered  in  three  concentric  ducts  corresponding  to  the 
center  core  igniter,  the  main  charge  and  the  external  ullage.  In  spite  of  the 
greater  completeness  of  modeling  of  the  configuration  of  the  bag  charge,  in 
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the  CALSPAN  code,  the  NOVA  code  has  been  shown  to  agree  better  with  a  stand¬ 
ard  155nun  center  core  configuration^^.  However,  neither  code  has  successfully 
been  able  to  predict  the  anomalous  behaviour  of  a  high  zone  base  ignited  155mm 
charge  except  by  ad  hoc  manipulation  of  the  data  base^^  or  by  assuming  grain 
fracture  due  to  impact  of  the  charge  on  the  base  of  the  pro;]  ectile^^  . 

In  order  to  determine  the  extent  to  which  the  quasi-one-dimensional  as¬ 
sumption  fundamental  to  the  NOVA  code  might  limit  its  predictive  capacity 
with  respect  to  bag  charges,  we  recently  performed  a  limited  theoretical 
study^^ .  By  modifying  the  code  to  accept  coupled  one  dimensional  models  of 
the  base  ignited  bag  charge  and  the  external  ullage  we  found  that  a  completely 
different  f lamespreading  path  could  be  produced  from  that  predicted  by  the  one 
dimensional  code  in  its  usual  form.  Primer  gas,  produced  in  the  breech,  can 
initiate  f lamespreading  in  the  rear  portion  of  the  charge  and,  by  flowing 
around  the  bag,  initiate  a  flame  at  the  front  so  that  the  charge  is  ignited 
by  two  flame  sheets. 

Thus  it  \^^as  concluded  that  predictions  of  f lamespreading  in  bag  charges 
may  be  seriously  in  error  unless  careful  attention  is  paid  to  the  precise 
distribution  of  all  ullage.  It  was  also  noted  that,  depending  on  its  permea¬ 
bility,  the  bag  material  may  play  an  extremely  influential  role  in  respect  to 
flamespreading.  Naturally,  the  history  of  the  longitudinal  pressure  wave  can¬ 
not  be  predicted  unless  the  flamaspreading  path  is  accurately  captured. 

We  may  summarize  the  situation  thus.  The  behaviour  and  amplitude  of 
longitudinal  pressure  waves  are  of  interest  because  of  their  correlation 
with  catastrophic  charge  malfunction  and  also  because  they  are  believed  to  be 
correlated  with  performance  variability^^.  Factors  known  to  influence  the 
wave  structure  are  the  igniter  venting  characteristics,  the  permeability  and 
progressivity  of  the  charge  and  the  detailed  distribution  of  axial  and  radial 
ullage.  It  is  also  thought  that  the  bag  material  and  liners  may  be  influen¬ 
tial.  One  dimensional  modeling  has  been  quite  successful  in  respect  to  the 
performance  of  case  charges  but  is  demonstrably  deficient  in  respect  to  bag 
charges , 


Accordingly,  the  prediction  of  pressure  wave  structure  in  bag  charges 
requires  at  least  a  quasi-two-dimensional  approach  with  careful  attention 
given  to  the  ullage  distribution  and  the  influence  of  the  bag  and  liner.  The 
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CALSPAN  model  does  admit  a  quasi-two-dimensional  formulation.  However,  the 
radial  expansion  of  the  bed  due  to  the  center  core  igniter  blast  is  not  con¬ 
sidered  and  the  influence  of  the  bag  is  neglected.  We  regard  these  as  serx- 
ous  physical  omissions.  In  addition,  the  numerical  analysis  contained  in  the 
CALSPAN  code  is  based  on  a  highly  diffusive,  and  therefore  inaccurate  scheme. 
Finally,  the  axial  distribution  of  ullage  is  represented  implicitly  and  a 
finite  difference  scheme  is  applied  in  spite  of  the  presence  of  discontinu¬ 
ities  in  the  porosity.  This  can  be  expected  to  lead  to  further  inaccuracy. 

Our  purpose  is  to  produce  an  axisymmetric  two  dimensional  code  which 
corrects  these  deficiencies  and  permits  a  significant  advance  in  our  ability 
to  predict  the  pressure  wave  amplitudes  associated  with  bag  charges.  It 
should  be  clearly  understood,  however,  that  while  our  investigation  into 
multi  dimensional  effects  will  be  a  major  undertaking,  it  will  not  close  all 
gaps  between  theory  and  experiment.  There  is  strong  evidence  that  many 
malfunctions^^ ’ may  involve  grain  fracture,  a  process  which  has  not  yet 
been  seriously  addressed  by  anyone  in  the  gun  modeling  community. 

The  present  report  addresses  the  formulation  of  the  model,  its  govern¬ 
ing  equations  and  the  method  of  solution.  Furthermore,  we  encode  the  model 
to  the  point  of  enabling  the  determination  of  numerical  solutions  for  a 
single  two  phase  flow  region  of  complex  two  dimensional  shape. 


1 . 3  Summary  of  Approach  and  Findings 

We  will  first  summarize  both  our  overall  approach  and  the  findings  of 
the  present  study.  Subsequently,  we  will  provide  a  more  detailed  discussion 
of  the  factors  which  have  influenced  the  approach  chosen  here  as  well  as  some 
of  the  implications  and  limitations  of  our  approach. 

Our  analytical  approach  may  be  summarized  as  follows.  We  divide  the  com¬ 
bustion  chamber  into  disjoint  regions  in  each  of  which  all  the  flow  variables 
may  be  regarded  as  continuously  differentiable.  In  particular,  one  such  re¬ 
gion  is  allocated  for  each  charge  increment  and  the  mixture  boundaries  always 
coincide  with  the  region  boundaries.  Accordingly,  a  precise  representation 
is  made  of  the  ullage.  The  flow  inhibition  associated  with  the  bag  material 
and  its  various  liners  may  be  embedded  accurately  as  boundary  conditions  link¬ 
ing  the  flow  in  one  region  with  that  in  its  neighbors.  For  flexibility  and 
economy  we  consider  that  the  flow  in  a  given  region  may  be  any  one  of  two 
dimensional  two  phase,  two  dimensional  single  phase,  quasi-one-dimensional  two 
phase,  quasi-one-dimensional  single  phase  or  lumped  parameter  single  phase. 

The  governing  equations  for  each  of  these  types  of  regions  consist  of 
balance  equations,  which  have  been  previously  derived!^) 26 ,  and  constitutive 
laws  for  which  it  is  necessary,  in  some  cases,  to  extrapolate  from  previous 
one  dimensional  laws.  The  approach  used  to  generate  numerical  solutions  is 
essentially  a  marching  technique  based  on  an  explicit  two  step  scheme.  All 
physical  regions  are  mapped,  in  a  time  dependent  manner,  onto  a  regular  com¬ 
putational  figure,  either  a  unit  line  or  a  unit  square. 
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The  numerical  approach  is  essentially  a  logical  extension  of  that  which 
we  have  previously  implemented  successfully  for  one  dimensional  problems^^. 

An  obvious  milestone  in  the  path  of  the  desired  extension  is  the  demonstra¬ 
tion  of  a  stable  integration  scheme  for  a  single  region  of  two  dimensional 
two  phase  flow  in  a  closed  chamber.  Moreover,  the  demonstrated  scheme  must 
be  clearly  capable  of  interfacing  in  a  natural  manner  with  the  additional 
requirements  imposed  by  the  considerations  of  multiple  regions,  including 
ullage,  separated  by  explicitly  represented  boundaries. 

The  specific  approach  that  we  take  to  this  particular  milestone  is  as 
follows.  We  accomodate  the  geometrical  complexity,  imposed  by  the  actual 
shape  of  a  typical  gun  propellant  chamber  and  an  intruding  projectile,  by 
mapping  the  physical  domain  onto  a  square.  An  equipotential  map  of  the  type 
used  by  Thompson  et  al^^  is  adopted  for  this  purpose.  The  solution  is  up¬ 
dated  at  interior  mesh  points  by  means  of  the  MacCormack  scheme^^  modified 
according  to  some  recommendations  of  Moretti^^.  The  boundary  values  are  up¬ 
dated  using  Kentzer’s^^  approach  to  the  method  of  characteristics,  also  modi¬ 
fied  in  accordance  with  some  recommendations  of  Moretti. 

We  present  stable  solutions  for  three  problems  based  on  nominal  data. 

The  first  two  differ  only  in  respect  to  the  representation  of  the  granular 
stress  term,  it  being  reversible  in  the  first  and  irreversible  in  the  second. 
The  container  is  taken  to  have  a  curved  breech  closure  similar  to  the  mushroom 
in  a  howitzer.  The  tube  is  tapered  and  the  base  of  the  projectile  is  flat. 
Subsequently,  in  the  third  problem  we  allow  the  projectile  to  have  a  boattail 
which  intrudes  into  the  two  phase  flow  producing  as  complicated  a  geometry  as 
we  expect  to  encounter^  for  any  single  region  of  continuous  two  phase  flow,  in 
future  applications  to  howitzer  charges. 

The  analytical  basis  for  the  model  as  a  whole  is  presented  in  section  2.0. 
In  section  3.0  we  present  the  details  of  the  method  of  solution  and  the  three 
solutions  described  in  the  preceding  paragraph. 

Having  now  summarized  our  approach  and  present  findings,  we  will  attempt 
to  clarify  matters  a  little  further  by  returning  to  the  physical  content  of 
the  model.  In  particular,  we  will  comment  on  the  nature  and  manner  of  model¬ 
ing  of  the  system  components  -  breech  block,  tube,  projectile,  charge  incre¬ 
ments,  igniter  increments  and  the  bag  cloth,  liner  and  additives. 
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Figures  1.2.1  and  1.2.2  have  respectively  illustrated  a  typical  case  and 
bag  charge  in  relation  to  the  tube,  breech  and  projectile.  The  points  that  we 
abstract  from  these  figures  are:  (a)  The  tube  is  tapered  and  has  fairly  sham 
corners  associated  with  the  forcing  cone  and  the  origin  of  rifling  where  the 
band  of  the  projectile  is  approximately  located  (b)  The  breech  of  the  Howitzer 
is  roughly  hemispherical,  (c)  The  boattail  of  the  projectile  has  a  rather  short 
intrusion  into  the  tube,  (d)  Rather  complex  packaging  elements  are  introduced 
at  the  front  end  of  the  case  charge,  (e)  Although  the  figure  does  not  make 
this  clear,  the  bag  charge  is  contained  in  a  rather  complex  package  consisting 
of  a  cloth  bag,  lined,  in  part,  with  lead  foil  to  assist  decoppering  of  the 
tube,  and  a  layer  of  talc  impregnated  cloth  to  assist  in  the  reduction  of 
erosion.  It  is  also  not  uncommon  for  a  bag  of  potassium  nitrate  to  be  stitched 
onto  one  or  both  ends  of  the  bag  in  order  to  eliminate  muzzle  flash,  (f)  It  is 
evident  that  the  bag  charge  is  not  an  axisymmetric  configuration.  Thus  our 
model  involves  a  geometrical  approximation  from  the  very  beginning.  While  it 
is  not  so  obvious,  the  same  is  true  of  the  case  charge  since  the  bayonet  primer 
vents  at  a  number  of  small  holes  and  thereby  introduces  a  three  dimensional 
flow. 


Each  region  occupied  by  propellant  is  modeled  as  a  heterogeneous  two  phase 
mixture  for  which  the  balance  equations  have  been  established^’  .  Thus,  as  in 
previous  work-^^  we  describe  the  flow  by  reference  to  macroscopic  state  var- 
iables  which  are  presumed  to  be  averages  formed  over  a  region  large  enough  to 
contain  many  particles  and  yet  small  by  comparison  with  the  overall  dimensions 
of  the  chamber.  It  should,  however,  be  noted  that  the  fundamental  assumption 
that  the  scale  of  heterogeneity  of  the  mixture  is  infinitesimal  by  comparison 
with  the  proportions  of  the  chamber  is  on  less  firm  ground  when  we  turn  from 
a  consideration  of  axial  to  radial  flow. 

Also,  as  in  previous  work,  the  balance  equations  consist  of  statements 
of  conservation  of  mass,  momentum  and  energy  for  the  gas  phase  and  of  mass  and 
momentum  for  the  solid  phase.  The  solid  phase  is  assumed  to  be  microscopically 
incompressible  so  that  an  energy  balance  is  tantamount  to  a  heat  transfer  re¬ 
lation.  Naturally,  the  bulk  average  temperature  of  the  solid  phase  is  not  of 
interest;  we  assume  that  ignition  is  predicated  on  the  surface  temperature. 

As  in  previous  work,  we  assume  that  drag,  heat  transfer  and  combustion  can  be 
represented  by  empirical  correlations  based  on  the  macroscopic  states  of  the 
two  phases.  The  motion  of  the  solid  phase  is  taken  to  depend  on  the  interphase 
drag,  the  macroscopic  pressure  gradient  in  the  gas  and  also  on  intergranular 
stresses.  These  latter  stresses  are  taken  to  depend  on  the  porosity  in  an  ir¬ 
reversible  manner  since  the  bed  modulus  is  known  to  be  much  greater  during  un¬ 
loading  than  during  loading.  In  neither  phase  do  we  consider  a  macroscopic 
shear  stress.  The  significant  effects  of  viscosity  and  of  heat  conduction 
are  assumed  to  be  confined  to  the  boundary  layers  around  the  particles  and  to 
be  embedded  in  the  correlations.  They  therefore  appear  in  the  balance  equa¬ 
tions  as  algebraic  non-homogeneous  terras  rather  than  as  higher  order  differ¬ 
entials. 

Each  region  occupied  by  the  mixture  is  treated  as  a  continuum  at  all 
stages  of  the  calculation.  However,  depending  principally  on  its  radial  ex¬ 
tent,  it  may  be  represented  as  either  fully  two  dimensional  or  quasi-one- 
dimensional.  Thus,  in  treating  the  bag  charge  shown  in  figure  1.2.2,  we  would 
expect  the  inner  mixture  region  -  the  igniter  -  to  be  treated  as  quasi-one- 
dimensional  at  most.  We  shall  have  further  comments  below,  on  the  modeling 
of  the  igniter. 
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In  addition,  it  is  expected  on  the  basis  of  studies  of  the  radial  flow 
that  following  the  completion  of  flamespreading,  no  mixture  region  will  re¬ 
quire  anything  more  detailed  than  a  quasi-one-dimensional  representation. 

The  only  non-trivial  radial  distribution,  once  flamespreading  is  complete,  is 
expected  to  be  that  of  porosity.  Thus  we  expect  that  even  if  a  fully  two  di¬ 
mensional  approach  is  adopted  during  flamespreading,  the  analysis  will  eventu¬ 
ally  be  amenable  to  a  quasi-two-dimensional  treatment  and,  possibly,  if  the 
radial  ullage  disappears  altogether,  the  latter  portions  of  the  interior  bal¬ 
listic  cycle  may  be  studied  by  the  conventional  quasi-one-dimensional  model. 

Similar  considerations  pertain  to  the  regions  which  we  employ  in  the  re¬ 
presentation  of  ullage.  These  are  all  presumed  to  contain  an  inviscid,  non¬ 
heat-conducting,  compressible  gas.  The  neglect  of  diffusion  is  in  accord  with 
the  purpose  of  the  model  which  is  to  predict  the  structure  of  longitudinal 
pressure  waves.  It  is  believed  that,  at  least  in  conventional  charges,  the 
influence  of  the  tube  mechanical  and  thermal  boundary  layers  will  be  negligible. 

Of  course,  if  we  were  interested  in  erosion  of  the  tube  the  behaviour  of 
these  boundary  layers  would  be  of  paramount  interest.  We  also  note  that  pre¬ 
dictions  of  maximum  pressure  of  the  tube  are  barely  influenced  by  the  heat 
loss  to  the  tube  and  even  the  muzzle  velocity  is  only  affected  to  the  extent 
of  3-5%.  There  is  no  doubt  in  our  mind  that  the  influence  of  wall  friction 
and  heat  loss  to  the  tube  by  the  gas  phase  will  be  negligible  in  so  far  as  the 
structure  of  the  longitudinal  pressure  waves  is  concerned.  An  additional  ef¬ 
fect  of  the  boundary  layer,  namely  flow  displacement  may  be  influential  in 
chambers  which  possess  significant  chambrage.  Separation  of  the  boundary 
layer  at  a  sharp  entrance  to  the  tube  could  conceivably  provide  an  appreciable 
constriction  of  the  inviscid  core  flow.  As  far  as  displacement  due  to  the 
transient  developing  boundary  layer  in  the  tube  is  concerned,  estimates  based 

Q  O 

on  the  momentum  integral  approach  of  Shelton  et  al  show  that  it  may  be  ne¬ 
glected  in  medium  caliber  weapons. 

The  regions  of  ullage  may  be  treated  as  two  dimensional,  quasi-one-dimen- 
sional  in  either  the  radial  or  the  axial  direction  or  as  lumped  parameter.  As 
with  the  mixture  regions,  the  principal  criterion  is  extent  in  the  radial  and 
axial  directions. 

The  influence  of  the  bag  material  is  readily  seen  to  be  embedded  in  the 
boundary  conditions  which  govern  the  mass  transfer  between  the  mixture  regions 
and  the  ullage  regions.  It  is  important  to  note  that  as  flamespreading  occurs 
the  mixture  regions  will  deform.  Our  approach  follows  this  deformation  pre¬ 
cisely  and  the  computational  regions  move  and  deform  so  that  the  boundaries 
are  always  located  precisely. 

We  now  summarize  the  representation  of  the  system  elements. 

Breech  Block,  Tube  and  Projectile 


The  configuration  of  the  chamber  as  defined  by  these  elements  is  repre¬ 
sented  accurately  through  tabular  input.  However,  corners  in  the  individual 
boundary  elements  will  not,  in  general,  coincide  with  continuum  mesh  points. 
__ 

'SheJXon,  5.,  A.  and  Saka,  P.  ”Study  o(^  Head:  TKan^^eA  and 
EAo^don  dn  Gun  EoAAelA^'  AfATL-TR-73-69  1973 
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The  projectile  is  assumed  to  move  as  a  rigid  body  and  to  be  opposed  by  a  pre¬ 
determined  law  of  barrel  resistance  due  to  engraving  of  the  rotating  band.  In 
general,  the  boundaries  of  the  chamber  are  impermeable  to  both  phases  but  gas 
loss  due  to  poor  obturation  in  the  breech  or  at  the  rotating  band  can  be  con¬ 
sidered.  Modeling  of  the  ignition  train  may  also  eventually  require  the  con¬ 
sideration  of  a  permeable  boundary  as  discussed  below. 

Propellant  Charge 


It  is  assumed  that  each  charge  element  consists  of  a  single  type  of  pro¬ 
pellant  and  that  within  each  such  element,  all  the  grains  are  initially  iden¬ 
tical,  However,  different  elements  need  not  contain  the  same  propellant.  The 
propellant  is  assumed  to  be  characterized  by  the  same  data  base  considered  in 
the  one  dimensional  model^®.  The  propellant  is  assumed  to  be  heated  by  the 
gas  phase  until  the  surface  temperature  exceeds  a  value  required  for  ignition 
and,  subsequently,  to  regress  in  accordance  with  a  steady  state  combustion  law. 
The  propellant  moves  as  a  consequence  of  drag,  pressure  gradient  and  also,  in¬ 
tergranular  stress. 

Propellant  Packaging 

Since  our  interest  is  primarily  directed  towards  bag  charges  we  do  not 
provide  a  model  of  the  closure  elements  typical  of  case  charges,  figure  1.2.1, 
and  which  have  played  an  important  role  in  our  one  dimensional  studieslO.  How¬ 
ever,  we  assume  that  each  charge  element  is  surrounded  by  a  bag  material  of 
negligible  thickness  and  mechanical  strength  and  which  has  variable  permeabil¬ 
ity  and  inertia.  The  permeability  is  assumed  to  depend  on  the  location  of  the 
surface  element  of  the  charge.  If  data  become  available,  the  permeability  may 
be  made  time  dependent  by  reference  to  the  local  mechanical  and  thermal  environ¬ 
ment.  Once  the  bag  material  is  determined  to  have  been  ruptured  or  rendered 
fully  permeable  it  disappears  from  consideration  in  the  model. 

Ignition  Train 

The  ignition  train  can  be  modeled  by  a  combination  of  an  external  source 
and  the  representation  of  the  center  core  charge  as  a  quasi-one-dimens ional 
two  phase  region.  However,  the  author  is  not  aware  of  data  which  demonstrate 
the  reliability  of  predictions  of  flamespreading  in  black  powder.  Therefore, 
the  representation  of  the  ignition  train  by  a  predetermined  rate  of  injection 
of  energetic  gas  is  thought  to  represent  a  more  fruitful  approach.  Some  com¬ 
parisons  of  theoretical  and  observed  flamespreading  in  black  powder  would  soon 
resolve  this  issue.  If  the  modeling  of  black  powder  combustion  is  indeed 
fruitful,  then  a  mixture  region  can,  in  principle,  be  allocated  for  the  re¬ 
presentation  of  a  black  powder  base  pad.  In  such  a  case,  the  breech  would  be 
taken  to  be  gas  permeable  in  order  to  permit  the  simulation  of  the  spit  hole. 
However,  that  level  of  modeling  is  not  recognized  by  the  present  study. 

As  regards  the  modeling  of  the  ignition  stimulus  by  an  external  source 
term,  we  note  that  previous  one  dimensional  studies  have  incorporated  a  cor¬ 
rection  for  the  volume  originally  occupied  by  the  condensed  phase  of  the 
igniter  ^*^5  1^5  26^  Here  we  take  the  attitude  that,  in  the  context  of  a  two 
dimensional  model,  the  correction  for  the  volume  of  the  condensed  phase  is  no 
longer  minor.  Therefore,  if  it  is  to  be  considered  at  all,  it  must  be  em¬ 
bedded  into  a  properly  posed  differential  equation  expressing  macroscopic  con¬ 
tinuity  of  the  substance  in  question. 
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We  conclude  the  discussion  of  the  model  with  some  comments  on  its  lirni- 
tations  of  applicability.  As  has  already  been  discussed,  the  neglect  of  the 
tube  boundary  layer  implies  that  the  model  cannot  be  applied  directly  to  prob¬ 
lems  of  tube  heat  transfer  and  erosion.  However,  the  model  may  be  used  to 
drive  a  boundary  layer  analysis. 

The  combustion  model  is  doubtless  oversimplified.  However,  the  present 
status  of  the  theory  of  transient  combustion  in  the  presence  of  cross  flow  is 
not  sufficiently  advanced  to  justify  a  more  complete  model.  Likewise,  the 
neglect  of  the  possibility  of  grain  fracture  may  well  prevent  the  model  from 
predicting  castastrophic  overpressures. 

We  also  note  that  the  scale  of  heterogeneity  of  medium  caliber  charges  is 
not  negligible  by  comparison  with  the  radial  dimensions  of  the  tube.  Accord¬ 
ingly,  attempts  to  resolve  too  many  of  the  details  of  the  flow  predicated  on 
the  ideal  assumptions  of  a  continuum  theory  are  likely  to  amount  to  self- 
deception.  Therefore,  the  liberal  use  of  quasi-one-dimensional  and  lumped 
parameter  representations  is  desirable  not  only  to  simplify  the  numerical  an¬ 
alysis  but  also  to  pose  the  simulation  consistently  with  the  underlying  assump¬ 
tions  . 

We  comment  now  on  the  choice  of  computational  examples  used  to  test  the 
numerical  scheme  for  two  dimensional  two  phase  flow.  Figure  1.2.2  shows  the 
bag  configuration  to  be  initially  rectangular  in  cross  section.  However,  we 
will  represent  the  charge  as  uniformly  occupying  the  entire  chamber,  includ¬ 
ing  the  region  around  the  boattail.  As  we  discuss  further  in  section  4.0, 
such  charge  configurations  are  of  immediate  interest  and  are  found  in  cased 
ammunition  for  which  the  projectile  intrusion  may  be  very  large.  The  rele¬ 
vance  of  such  a  computational  example  to  problems  involving  bag  charges  is 
simply  this:  as  a  consequence  of  igniter  blast,  the  bed  will  expand  and  the 
most  complicated  geometry  we  expect  ,  for  any  single  region  of  continuous  flow, 
corresponds  to  that  defined  by  the  external  boundaries. 

We  conclude  our  introduction  with  some  comments  on  the  choice  of  the 
method  of  solution. 

In  terms  of  the  gun  problem  and  convective  f ],amespreading  in  general,  we 
note  that  several  approaches  have  been  used  with  reasonable  success.  We  have 
used  a  variation  of  the  MacCormack  scheme^^  in  our  most  recent  work^^  supple¬ 
mented  by  the  method  of  characteristics  at  the  external  and  internal  boun¬ 
daries.  Kuo  has  used  a  Lax-Wendroff  type  scheme  also  supplemented  by  the 
method  of  characteristics^^ .  Krier  has  used  the  Richtmyer  scheme  with  arti¬ 
ficial  viscosity  and  Fisher  has  used  the  Lax  scheme^ ^ ^  Xt  should  be  noted 
that  the  Lax  scheme  is  very  inaccurate  due  to  the  large  implicit  viscosity  em¬ 
bedded  in  the  differencing  technique.  Thus  stability  is  acquired  only  with  a 
loss  of  accuracy.  The  artificial  viscosity  used  by  Krier  and  to  a  lesser  ex¬ 
tent  by  Kuo  does  not  appear  to  be  necessary  in  problems  of  gun  flamespreading 
according  to  our  own  experience.  Because  we  have  bad  satisfactory  results 
using  the  MacCormack  scheme  and  the  method  of  characteristics  in  our  one  di¬ 
mensional  simulation,  it  is  natural  for  us  to  have  considered  using  the  same 
approach  in  our  two  dimensional  simulations. 


However,  many  papers  are  available  to  describe  solutions  of  two  dimen¬ 
sional  unsteady  multiphase  flow,  almost  all  being  the  work  of  the  Los  Alam.os 
group.  Fortunately,  there  appears  to  be  complete  compatibility  between  our 
balance  equations  and  theirs so  that  their  methods  are  of  particular  in¬ 
terest  . 

We  have  previously  commented  on  two  papers  by  Harlow  and  Amsden^^’^^. 

These  papers  demonstrate  the  feasibility  of  obtaining  numerical  solutions  in 
two  dimensions  with  several  internal  boundaries,  all  represented  implicitly. 
Several  programs  are  documented  by  the  Los  Alamos  group  including  K-TIF^^, 
K-FIX^S  and  KACHINA^^. 

The  methodology  represented  by  these  papers  suffers  from  three  difficul¬ 
ties  when  we  relate  it  to  our  own  application.  In  the  first  place,  the  im¬ 
plicit  finite  difference  scheme  is  rather  complicated.  This  is  not  a  serious 
objection  per  se.  However,  if  modifications  are  required  in  our  application, 
the  numerical  ramifications  may  well  be  severe.  For  example,  Lee  et  al^*^  re¬ 
port  difficulties  experienced  with  the  KACHINA  code  in  applications  to  reac¬ 
tor  containment  problems.  These  difficulties  included  loss  of  stability  and, 
not  surprisingly,  severe  numerical  diffusion  problems.  Secondly,  there  is  a 
surprising  lack  of  regard  even  for  external  boundary  conditions  in  these  papers. 
As  usual,  the  reflection  technique,  which  has  been  criticized  by  Moretti^f  is 
used  to  describe  impermeable  boundaries.  More  seriously,  a  simple  continua- 
tive  condition  is  used  in  KACHINA  to  describe  outflow,  even  for  a  subsonic  con¬ 
dition.  Since  the  outflow  depends  upon  the  external  pressure  in  the  subsonic 
mode,  it  is  difficult  to  know  what  such  solutions  mean.  Finally,  the  failure 
to  capture  explicitly  the  internal  boundaries  defined  by  jumps  in  porosity  is 

a  very  serious  omission  as  far  as  our  application  is  concerned. 
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The  boundary  conditions  at  the  internal  surfaces  drive  the  flamespread¬ 
ing  to  a  large  degree  since  it  is  through  these  conditions  that  we  express 
the  behaviour  of  the  bag  and  liner  materials.  Moreover,  even  a  very  small 
gap  around  the  charge  can  be  extremely  influential  in  respect  to  flamespread¬ 
ing  and  the  subsequent  structure  of  the  longitudinal  pressure  waves^^.  Accord¬ 
ingly,  the  prospect  of  numerical  diffusion  which  would  obscure  the  existence  of 
ullage  is  inadmissible  in  our  application. 

We  close  with  the  following  general  observations  on  the.  relative  merits  of 
explicit  and  implicit  schemes.  Explicit  schemes  are  attractive  because  of  their 
computational  simplicity.  Implicit  schemes  are  generally  more  complex,  requir¬ 
ing  greater  computation  per  time  step,  but  offer  stability  over  time  steps  which 
may  greatly  exceed  those  allowable  with  explicit  schemes.  When  the  allowable 
time  step  is  sufficiently  greater  than  that  obtainable  with  an  explicit  scheme, 
the  additional  computational  burden  per  time  step,  of  the  implicit  scheme,  be¬ 
comes  economically  viable.  Indeed,  for  problems  which  involve  diffusion,  the 
implicit  schemes  offer  significant  overall  economy  relative  to  explicit  schemes. 

In  the  present  case,  however,  the  nonhomogeneous  terms  are  known  to  exert 
a  dominant  influence  on  the  time  step,  at  least  during  flamespreading.  This  has 
been^not  only  our  own  experience  using  an  explicit  scheme^ but  also  that  of 
Kuo^  who  used  an  implicit  scheme  in  a  pioneering  study  of  convective  flamespread¬ 
ing  in  a  stationary  porous  bed.  Accordingly,  as  the  explicit  and  implicit 
approaches  to  a  finite  difference  scheme  are  expected  to  be  similarly  constrained, 
in  regard  to  time  step,  the  computational  simplicity  of  the  explicit  scheme 
weighs  strongly  in  its  favor. 
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2.0  GOVERNING  EQUATIONS 


As  we  have  discussed  in  the  introduction,  our  approach  to  the  analysis 
of  flamespreading  and  pressure  wave  development  in  a  complex  propelling 
charge  is  based  on  the  solution  of  the  balance  equations  for  each  region  of 
continuous  flow  together  with  explicit  boundary  conditions  which  reflect  the 
couplings  between  the  regions.  Our  purpose,  in  the  present  section,  is  to 
take  note  of  these  equations  in  sufficient  detail  as  to  provide  an  analyti¬ 
cal  basis  for  our  model.  At  the  same  time,  we  will  define  that  subset  of  the 
governing  equations  whose  numerical  solution  is  undertaken  in  the  present 
study. 

Thus,  section  2.1  includes  equations  for  both  two  phase  and  single  phase 
flow  in  both  one  and  two  dimensional  regions.  However,  the  subsequent  numer¬ 
ical  studies  pertain  solely  to  a  two  dimensional  region  of  two  phase  flow. 

In  section  2.2  we  discuss  the  constitutive  laws  required  to  close  the  balance 
equations.  In  general,  the  extension  from  earlier  one  dimensional  studies  is 
analytically  straightforward.  However,  new  restrictions  on  the  validity  of 
the  correlations  for  the  interphase  processes  arise  and  are  duly  noted.  We 
also  note  the  simplified  forms  of  the  constitutive  laws  used  in  the  present 
numerical  studies . 

Following  the  discussion  of  the  constitutive  laws,  we  consider,  in  sec¬ 
tion  2.3,  the  conditions  which  apply  at  both  the  external  boundaries  of  the 
computauional  domain  and  at  the  internal  boundaries  which  separate  regions 
of  continuous  flow.  As  part  of  this  topic  we  discuss  the  behaviour  of  the 
bag  or  packaging  material  and  we  note  the  form  of  the  constitutive  laws  suit¬ 
able  for  the  description  of  bag  strength,  inertia  and  permeability.  We  also 
comment  on  the  initial  conditions. 

We  conclude,  in  section  2.4,  with  a  characteristic  analysis  of  the  bal¬ 
ance  equations.  The  resulting  conditions  of  compatibility  play  an  essential 
role  in  respect  to  the  method  of  solution,  as  we  discuss  further  in  section 
3.0.  Results  are  developed  for  both  single  and  two  phase  two  dimensional 
flow,  although  only  the  latter  are  used  in  the  present  numerical  studies.  We 
also  note  that  the  discussion  of  section  2.4  anticipates  the  numerical  tech¬ 
nique  and  therefore  accounts  for  a  transformation  of  independent  variables 
whose  purpose  is  to  map  each  region  of  continuous  flow  onto  a  regular  figure. 
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2. 1  Systems  of  Balance  Equations 


We  require  statements  of  the  balance  equations  for  both  two  dimensional 
and  quasi-one-dimensional  flow  of  either  a  single  or  two  phase  substance.  In 
the  latter  case  we  formulate  the  balance  equations  for  each  of  the  individual 
species.  We  also  consider  a  lumped  parameter  description  for  a  region  of 
single  phase  flow,  it  being  understood  throughout  this  report  that  the  single 
phase  flow  is  always  that  defined  by  the  gas  phase. 

The  two  phase  flow  balance  equations  are  as  given  previously  and  are 
governing  equations  for  macroscopic  properties  of  each  of  the  phases  or  aver¬ 
ages  formed  over  regions  large  in  comparison  with  the  scale  of  heterogeneity 
of  the  mixture.  As  in  our  previous  work,  the  complex  boundary  layer  phenomena - 
drag,  heat  transfer,  combustion  rate  -  appear  as  nonhomogeneous  terms  or  alge¬ 
braic  entities  in  the  balance  equations.  It  is  assumed  that  empirical  corre¬ 
lations  are  available  to  relate  these  processes  to  the  macroscopic  state 
variables  as  we  describe  further  in  section  2.2. 

The  quasi-one-dimens ional  formulation  of  the  two  phase  flow  equations 
contains  a  provision  for  mass  exchange  with  neighboring  regions.  Both  the 
quasi-one-dimensional  and  the  two  dimensional  systems  incorporate  a  source 
term  to  reflect  an  ignition  stimulus.  However,  heat  loss  to  the  tube  is  ne¬ 
glected  and  the  stress  tensor  is  taken  to  be  isotropic  for  both  species  so 
that  resistance  to  shearing  is  not  considered  for  either  the  gas  or  the  solid 
phase. 

We  assume  the  single  phase  continuum  flow  to  be  inviscid  and  non-heat 
conducting.  However,  v/e  do  retain  the  source  terms  which  may  embed  either  an 
ignition  stimulus  or  mass  transfer  from  a  neighboring  region, 

2.1.1  Two  Dimensional  Two  Phase  Flow 


In  cylindrical  coordinates  such  that  z  is  the  axial  coordinate ,  r  is  the 
radial  coordinate  and  t  is  the  time,  the  balance  equations  take  the  forms: 

Balance  of  Mass  of  Gas  Phase 


^  . 

Dt 


r3u 

+  37I 


=  m 


+  i|,. 


cpv 

r 


2.1.1. 1 


The  notation  conforms  with  that  used  previously.  We  have  p,  the  density  of 
the  gas,  £  the  porosity,  u  and  v  the  z-  and  r-  components  of  gas  phase  vel¬ 
ocity,  D/Dt  the  convective  derivative  along  the  gas  phase  streamline,  4;  the 
source  term  associated  with  a  stimulus,  m  the  rate  of  production  of  gas  due 
to  combustion  of  the  solid  phase.  It  should  be  noted,  in  accordance  with  the 
discussion  given  in  the  introduction,  that  we  neglect  the  volume  occupied  by 
the  condensed  phase  of  the  ignition  stimulus. 

We  recall: 


D 

Dt 


9 

9t 


+  U  +  V 

9z 


9r 


2. 1.1.2 
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2. 1.1.3 


m  =  (1-e)  dp  =  s  dp 
V  p  P  P 

Here  Sp,  Vp  are  the  surface  area  and  volume  of  an  individual  grain  and  d  is 
the  rate  of  surface  regression.  We  have  introduced  Sp  as  the  surface  area 
per  unit  volume. 

Balance  of  Momentum  of  Gas  Phase 


Du  ,  „ 


-  f  +  m(u  -  u)  -  iIju 
P 


2. 1.1.4 


Here  u  is  the  velocity  with  ^pomponents  u  and  v.  Also,  in  cylindrical  coord¬ 
inates  Vp  =  (9p/8z,  9p/9r).  The  vector  form  of  the  interphase  drag  f  should 
be  noted. 

Balance  of  Energy  of  Gas  Phase 


^  +  ^]  +  P  ^  -  f  • 


Dt 


P 


u-u 

•  p  I  L 

-  s  q  +  m  (e  -e  +  — ■ 

P  P  Pp  2g^ 


.  J  /  1  U*U.  V 

^  -  "P  r 


2. 1.1.5 


Here  e  =  e(p,p)  is  the  internal  energy  of  the  gas  phase  and  q  is  the  inter¬ 
phase  heat  transfer  per  unit  surface  area  of  the  solid  phase. 

Balance  of  Mass  of  Solid  Phase 


De 


Dt 

P 


(l-£)[ 


9u 

9z 


+ 


V 

(l-e) 

r 


2. 1.1.6 


The  subscript  p  denotes  properties  of  the  solid  phase  and  D/Dt  is  defined 
by  analogy  with  2. 1.1. 2.  ^ 

Balance  of  Momentum  of  Solid  Phase 

Du 

(l-e)p  5^  +  (l-e)g^Vp  +  g^Va  =  f  2. 1.1.7 

P  p 

The  vector  form  of  this  equation  should  be  noted.  We  have  o'  =  (1-£)R(£,£) 
where  R  is  the  average  stress  due  to  contacts  between  particles  and 
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will  be  assumed  to  depend  upon  porosity  in  an  irreversible  manner. 
2.1.2  Quasi-One-Dimensional  Two  Phase  Flow 


We  assume  that  for  the  applications  of  interest  to  us,  the  non- trivial 
spacewise  coordinate  is  aligned  with  the  axis  of  the  tube  and  that  u,  Up  are 
the  non-trivial  components  of  gas  and  solid  phase  velocity  respectively.  The 
cross  sectional  area  of  the  annulus  through  which  the  flow  occurs  is  taken  to 
be  A(2,t)  and  therefore  depends  upon  both  position  and  time.  It  is  supposed 
that  the  circumferential  boundaries  are  permeable  to  the  gas  phase  and  that 
mass  transfers  must  be  considered.  We  use  and  to  denote  respectively 
the  radii  of  a  circumferential  surface  on  which  influx  (mj^)  or  efflux  (mo) 
occur.  Attention  should  be  paid  to  this  convention.  The  subscripts  i  and  o 
do  not  refer  to  the  interior  and  exterior  surfaces,  only  to  the  direction  of 
mass  transfer.  We  understand  m^  and  m^  to  represent  rates  of  transfer  per 
unit  surface  area.  Moreover,  we  will  also  denote  the  properties  transported 
with  m^  by  the  subscript  i.  Thus  u^  will  be  the  axial  velocity  associated 
with  the  incoming  gas.  The  exiting  properties  are,  of  course,  those  of  the 
gas  in  the  quasi-one-dimensional  region  presently  under  consideration. 

Balance  of  Mass  of  Gas  Phase 


r\  Pi 

^  sAp  -H  —  eApu  =  Am  +  AiIj  +  27r[ZR.m.  ™  ZR  m  ]  2.  1.2.1 

dt  ^  dz  ^  ^  ^11  o  o-* 

Of  course,  the  summations  are  over  all  entering  and  all  exiting  fluxes - 
Balance  of  Momentum  of  Gas  Phase 
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£P  +  eg 


O  dz 


f  -  l|ju  H-  m(u  -u)  + 
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.  ZR.m. (u.-u) 

A  111 


2. 1.2.2 


Balance  of  Energy  of  Gas  Phase 


De  p  DeA  ,  9u  f  .  . 


r  P 

+  m  e  -  e  +  i 


(u-u„)2 


P  Pp 
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2^  P  V*  D 
—  ™  Zm  R 
A  p  o  o 


2. 1.2.3 
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Balance  of  Mass  of  Solid  Phase 


f  5^  (1-£)A+  (l-e)  ^ 

p 


2, 1.2.4 


Balance  of  Momentum  of  Solid  Phase 


%  ^  Dt  9z  9z 

P 


2. 1.2.5 


2.1.3  Two  Dimensional  Single  Phase  (Gas)  Flow 

These  equations  are  quite  familiar.  We  also  note  that  they  represent 
the  limiting  forms  of  2. 1.1.1,  2. 1.1.4  and  2. 1.1,5  as  £->^1,  bearing  in  mind 
that  f,  m  and  q->o.  We  have: 

Balance  of  Mass 


Dp  .  9^1  -  I  pv 

m  -  r 


2. 1.3.1 


Balance  of  Momentum 


Du 

P  ^  +  g^Vp 


2. 1.3.2 


Balance  of  Enerc 


-9u 

■9z  9rJ 


,  r  ,  U*U  V 

■  JT  ■  '1  ■  ■■  ? 


2. 1.3.3 


2.1.4  Quasi-Qne~Dimensional  Single  Phase  (Gas)  Flow 

As  in  the  previous  section  the  balance  equations  for  this  case  follow 
from  those  set  forth  in  section  2.1.2  by  taking  the  limit  e  0.  We  have, 
assuming  the  axial  direction  to  be  non-trival: 

Balance  of  Mass 


7^  Ap  +  7^  Apu  =  Aij;  +  2Tr[ZR.m.-  ZR  m  ] 
9t  ^  9z  ^  ^  *-11  o  0-* 


2. 1.4.1 
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Balance  of  Momentum 


^  Du  ^ 

P  57  +  «o  sf 
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2. 1.4.2 


Balance  of  Energy 
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2. 1.4.3 


If,  in  fact,  the  non-trivial  direction  is  not  axial,  the  divergence  and 
transfer  terras  raust  be  suitably  modified. 

2.1.5  Lumped  Parameter  Single  Phase  (Gas)  Flow 

As  in  previous  work  we  provide  balance  equations  only  for  mass  and 
energy.  It  is  assumed  that  the  velocity  of  the  gas  in  the  lumped  parameter 
region  can  be  deduced  from  the  boundary  values  by  interpolation.  Using  V 
to  denote  the  volume  of  the  region  and  S  the  bounding  surface  we  have: 


dV  j 

- —  =  da 

dt 


2. 1.5. 1 


where  w  is  the  boundary  velocity  and  n  is  the  outward  facing  normal.  The 
mass  balance  is : 


Jip  dv  +  Em. 
V  ^ 


Em 

o 


2. 1.5.2 


•  • 

where  the  m-  and  m^  now  refer  to  the  total  fluxes  rather  than  the  fluxes  per 
10  ^ 
unit  area  used  previously. 
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The  energy  balance  is: 


where 


dt 


pEV 


dv  +  Zm. 

lb  1 


—  Em  (E  +  ■^) 

o  p' 


is  the  total  energy  e  +  vp-jlg 
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2.2  Constitutive  Laws 


The  constitutive  laws  required  to  close  the  foregoing  systems  of  bal¬ 
ance  equations  are  of  several  types.  As  is  usual  in  fluid  mechanical  prob¬ 
lems  we  require  an  equation  of  state  for  the  gas.  However,  we  also  need 
formulae  for  the  particle  surface  area  and  volume  and  for  the  granular  stress 
due  to  contacts  among  particles.  In  addition  to  these,  we  require  formulae 
for  the  interphase  transfer  processes  -  drag,  heat  transfer  and  combustion. 

In  order  to  discuss  combustion  we  will  have  to  introduce  the  surface  temper¬ 
ature  of  the  solid  phase  and  a  criterion  for  ignition.  All  of  the  above 
topics  are  addressed  in  the  present  section.  However,  the  mass  transfers 
from  one  region  to  another,  due  to  the  gas  permeability  of  the  boundary  be¬ 
tween  them  will  be  discussed  in  section  2.3. 

2.2.1  Equation  of  State  of  Gas 


It  is  assumed  that  the  gas  obeys  the  covolume  equation  of  state: 

e  =  c  T  =  2.2. l.l 

V  (y-1)p 

where  b  is  the  covolume,  y  is  the  ratio  of  specific  heats  and  c,,  is  the 
specific  heat  at  constant  volume.  In  previous  studies  we  have  considered 
the  composition  dependence  of  the  molecular  weight  and  the  specific  heats. 
From  an  analytical  viewpoint  this  simply  amounts  to  the  addition  of  two 
first  order  partial  differential  equations  to  the  model.  Physically,  the 
consequences  of  the  composition  dependence  are  found  to  be  small  as  the 
values  of  molecular  weight  and  specific  heat  appropriate  to  the  propellant 
are  quickly  attained.  Moreover,  we  do  not  consider,  herein,  mixtures  of 
propellants.  Thus,  in  the  present  work  we  treat  the  molecular  weight  and 
the  ratio  of  specific  heats  as  constants,  an  assumption  which  can  be  relaxed 
easily  when  future  studies  so  demand. 

2.2.2  Granular  Stress  Law 


The  formulation  of  the  granular  stress  law  is  of  considerable  importance 
not  only  because  of  the  physical  role  played  by  the  forces  transmitted  from 
grain  to  grain  but  also  because  of  the  influence  of  the  stress  law  on  the 
degree  of  hyperbolicity  of  the  balance  equations.  Thus,  although  the  present 
study  can  tolerate  simplifications  in  certain  of  the  constitutive  laws  for 
the  sake  of  temporary  expedience,  it  is  essential  that  we  evaluate  our  numer¬ 
ical  scheme  in  the  context  of  the  physically  complete  granular  stress  law. 


The  granular  stress  is  taken  to^de^end  on  porosity  and  also  on  the  di¬ 


rection  of  loading.  As  in  the  past 


we  embed  the  constitutive  law  into 


the  formula  for  the  rate  of  propagation  of  intergranular  disturbances: 


a(e) 


da-i'^ 

p  d£-^ 
p 


2. 2. 2.1 
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We  may  recast  2.2.2. 1  into  a  form  more  suitable  for  numerical  integra¬ 
tion,  namely; 


Da 

Dt 


P 


a^  D£ 


2. 2. 2. 2 


In  order  to  formulate  the  functional  behaviour  of  aCe)  we  introduce  Gq*  ^he 
settling  porosity  of  the  bed,  and  values  of  a(e)  equal  to  a^  and  3^2  which 
respectively  correspond  to  loading  at  Cq  and  to  unloading/reloading.  The 
nominal  loading  curve,  corresponding  to  monotonic  compaction  of  the  bed 
from  to  a  smaller  value  of  the  porosity  e  is  given  by: 


a 


a 

nom 


The  functional  dependence  of  a(e)  may  now  be  stated  as: 


2.2.2. 3 


na.e/e  ife<0,a  =  a  ,£<e 
1  o  —  "  nom  — 


a(£)  -  A 


a.  ifO<a<a  ,e<e 

2  ^  nom  —  o 

or  if  £  >  0,  a  =  a  ,  £  <  £ 
nom  —  o 


0 


if  a  =  0  and  £  >  0  or  if  e  >  £ 


2.2. 2.4 


Here  we  understand  £  to  mean  Ds/Dtp.  In  general,  it  appears  that  32  ^  3a ^ 
so  that  the  unload/reload  modulus  is  roughly  one  order  of  magnitude  greater 
than  the  nominal  loading  value.  We  emphasize  that  this  constitutive  law 
does  not  admit  any  granular  stress  when  £  >  £q,  that  is,  when  the  bed  is 
dispersed.  We  also  emphasize  the  absence  of  the  stress  decay  factor  used 
in  previous  work  for  purely  numerical  reasonsl4-18 ^ 

2.2.3  Propellant  Form  Functions 


It  is  assumed,  in  the  present  study,  that  the  particles  are  multi- 
perforated  cylinders  having  initial  length  Lq,  external  diameter  Dq  and  per¬ 
foration  diameter  dQ.  The  surface  area  per  unit  volume  is  related  to  the 
individual  surface  area  Sp  and  volume  Vp  of  the  particle  according  to: 

Sp  =  (l-£)Sp/Vp  2.2.3. 1 


Until  such  time  as  slivering  occurs,  that  is  to  say  the  time  at  which  the 
regressing  perforation  surfaces  intersect,  the  surface  area  and  volume 
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are  given  by: 


S  =  tt(L  -  2d)  [  (D  -  2d)  +  N(d  +  2d)] 

poo  o 

+  -it/2[(D  -  2d)2  -  N(d  +  2d)^]  2. 2. 3. 2 

o  o 

V  =  -rrCL  -  2d)[(D  -2d)^  -  N(d  +  2d)^]  2. 2. 3. 3 

poo  o 


where  N  is  the  number  of  perforations  and  d  is  the  total  linear  surface 
regression,  assumed  uniform  over  the  surface. 

Formulae  for  other  grain  geometries  are  simple  to  determine  pro¬ 
vided  that  the  initial  configuration  is  regular  and  that  slivering  has 
not  occurred.  Formulae  for  the  form  functions  following  the  slivering  of 
seven  perforation  grains  may  be  found  in  Krier  et  al^^ . 

2.2.4  Interphase  Drag 
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In  previous  work  we  have  assumed  the  interphase  drag  to  be  composed 
of  a  steady  state  component  and  a  transient  component  reflecting  the  virtual 
mass  effect: 


^  ^  Du 

f  -  -  jJij  2.2.4 

P 

where  3  is  the  virtual  mass  constant.  However,  while  our  one  dimensional 
modeling  exercises  have  carried  the  influence  of  virtual  mass  throughout 
the  analysis  and  coding,  actual  calculations  have  almost  always  been  based 
on  the  value  3  ”  0,  due  to  the  absence  of  reliable  independent  data  for 
the  actual  value.  This  being  the  case,  we  intend  to  neglect  the  virtual 
mass  effect  ab  initio  as  a  simplification  of  the  present  study. 

The  formula  for  the  steady  state  component  of  the  drag  is  assumed  to 
follow  as  the  obvious  generalization  of  the  one  dimensional  form: 
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where  Dp  is  t^he  effective  diameter  of  the  particles,  given  as  Dp  =  6Vp/Sp 
and  Rep  =  p|u-Up|Dp/|j  is  the  Reynolds  number.  In  previous  studies  we  have 
taken  Cpg  equal  to  the  value  reported  by  Ergun^^,  for  packed  beds  of  various 
shaped  particles,  modified  according  to  the  tortuosity  factor  of  Anderssen^^ 
as  the  bed  becomes  dispersed.  In  the  present  study  we  treat  Cpg  as  constant 
and  equal  to  Ergun’s  limiting  value  of  1,75  for  high  Reynolds  number  flows. 
Success  with  a  vector  form  of  Ergun’s  correlation  has  been  reported  by  Stanek 
and  Szekely^^. 


The  value  of  Cpg  =  1.75  for  a  packed  bed  was  determined  by  Ergun  for  a 
Reynolds  number  regime  which  is,  in  fact,  low  by  comparison  with  the  values 
expected  in  interior  ballistic  flows.  Kuo^^  has  shown  that  continues 
to  decrease  at  least  until  Rep  =  40000,  in  beds  of  ball  propellant  and  that 
a  value  of  Cpg  -'1. 3  may  be  more  typical  of  the  interior  ballistic  situation 
in  which  Rep  -10^.  In  computing  the  effective  diameter  of  the  particles, 
the  question  arises  as  to  the  actual  contribution  of  the  surface  area  of  the 
perforations.  Robbins^®  has  presented  data  which  indicate  that  Cpp  *^1.2  if 
the  entire  surface  area  is  used  to  compute  the  effective  diameter  and  that 
Cpp  ^1.9  if  only  the  contribution  of  the  external  surface  is  considered. 


We  also  note  that  relatively  few  grains  are  present  as  we  traverse  the 
tube  in  a  radial  direction.  Generally  speaking,  a  medium  caliber  weapon 
whose  diameter  is,  say,  15-20  cms  has  a  charge  consisting  of  grains  whose 
lengths  are  typically  2-3  cms  and  whose  diameters  are  typically  1  cm.  The 
reported  values  of  Cpp  are  based  on  flows  through  long  columns  for  which  en¬ 
trance  and  exit  effects  are  negligible.  Evidently,  as  few  as  5  grains  may 
lie  between  the  centerline  of  the  tube  and  the  inside  of  the  wall.  By  in¬ 
vestigating  the  drag  on  individual  particles,  van  der  Merwe  and  Gauvin^^ 
showed  that  seven  or  eight  layers  may  be  needed  to  produce  conditions  typical 
of  a  long  bed. 
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2.2.5  Interphase  Heat  Transfer 


50 

Previously  we  have  used  both  the  Denton  correlation  for  packed  beds 
and  the  Gelperin-Einstein^^  correlation  for  dispersed  and  packed  beds.  How¬ 
ever,  the  two  are  sufficiently  close  that  we  consider,  herein,  only  the  cor¬ 
relation  of  Gelperin-Einstein.  We  express  the  heat  transfer  in  the  form; 


Nu  =  0.4  2.2.5. 1 

P  P 


where  Nu  =  liD  /k  _ 

P  P  t 


Re  =  p.|u-u  Id  /jj 
p  P  P  f 


h  =  q/(T-Tp) 

The  subscript  f  denotes  an  evaluation  of  properties  at  the  film  temperature 
(T+Tp)/2  where  T  and  Tp  are  respectively  the  gas  bulk  average  temperature 
and  the  particle  surface  average  temperature.  The  viscosity  is  taken  to 
have  a  Sutherland  type  dependence  on  temperature; 


y 


0.134064 


(T/298)^'^ 
T  +  110 


2. 2. 5. 2 


The  thermal  conductivity  follows  from  the  Prandtl  number  which  is  assumed 
to  satisfy: 


k  9y-5 


2. 2.5. 3 


2.2.6  Solid  Phase  Surface  Temperature 


In  previous  one  dimensional  studies  we  have  treated  the  determination 
of  surface  temperature  as  though  conditions  on  the  surface  were  essentially 
homogeneous.  That  is  to  say,  the  heat  flux  given  by  2.2.5. 1  was  treated  as 
uniform  over  the  particle.  Moreover,  conduction  to  the  interior  was  assumed 
to  be  governed  by  a  one  dimensional  form  of  the  heat  diffusion  equation  and 
solved  by  means  of  a  cubic  profile^^ . 

*  Vmton,  W.  H.  "GeneAri/  V^cu^^yion  on  He.cU 
Jn6t.  knah.  Eng.  and  Am.  Soc.  Eng.  795/ 
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The  basis  for  the  assumption  of  uniformity  derives  from  observations  of 
quenched  beds .  In  these  the  surfaces  show  regular  regression  of  all  surfaces. 
Thus  it  is  inferred  that  surface  relaxation  phenomena  such  as  transverse 
flamespreading  proceed  very  quickly.  This  need  not  always  be  true,  however. 
Data  acquired  by  Jakus^^  in  conditions  of  low  pressure  ignition  revealed  very 
poor  uniformity.  Also  it  is  thought  that  the  assumption  of  surface  uniform¬ 
ity  may  be  poor  for  long  grains,  particularly  stick  or  cord,  which  have 
perforations . 

In  the  multi  dimensional  calculation  we  have,  in  principle,  an  essen¬ 
tially  new  feature.  The  stagnation  point  on  a  given  particle  may  vary  as  the 
direction  of  relative  flow  changes.  In  keeping  with  the  previous  assumption 
of  surface  uniformity,  this  feature  is  neglected.  Thus  the  surface  tempera¬ 
ture  is  given  by: 


T 

P 


hH 
k  2 
P 


[(I  - 

■o 


2  hH  ,  4  hTH 

3  3  ■ 

P  P 


where  T 


P 


o 


is 


the  initial  surface  temperature  and  H  satisfies: 


2,2.6, 1 


DH 

Dt 

P 


a  q 
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2.2. 6.2 


2,2.7  Ignition  and  Combustion 

Ignition  is  assumed  to  occur  when  the  surface  temperature  exceeds  a 
predetermined  value.  The  rate  of  surface  regression  is  given  by: 


Dd 


Dt 

P 


B  +  B 
1  2 


1.1.1 . 1 


It  should  be  noted  that  only  one  of  2. 2. 6. 2  and  2.2,7  A  has  to  be  solved 
at  each  point  according  as  the  temperature  is  less  than  or  equal  to  the 
ignition  temperature. 


2.2.8  Primer  Stimulus 


The  primer  stimulus  is  assumed  to  be  given  in  tabular  form  as  a 
predetermined  function  of  space  and  time. 
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2.3  Initial  and  Boundary  Conditions 


The  initial  conditions  are  more  or  less  self  evident  and  require  only 
minimal  attention.  For  problems  of  interest  in  interior  ballistics  we  sup¬ 
pose  that  both  phases  are  at  rest  and  at  atmospheric  pressure.  The  tempera¬ 
tures  of  the  two  phases  may  differ  but  are  uniform  throughout  each  of  the 
respective  media.  The  porosity  is  piecewise  continuous,  the  discontinuities 
being  defined  by  the  presence  of  ullage  or  of  the  boundary  between  bags. 

We  therefore  turn  our  subsequent  attention  to  the  analysis  of  the  bound¬ 
ary  conditions.  The  natural  starting  point  for  this  discussion  is  the  form 
of  the  jump  conditions  at  a  discontinuity  in  two  phase  flow.  This  topic  is 
addressed  in  section  2.3.1.  Then,  in  sections  2.3.2,  2.3.3  and  2.3.4  we  suc¬ 
cessively  consider  the  boundary  conditions  which  apply  at  the  external  boun¬ 
daries,  the  gas  permeable  internal  boundaries  and  at  the  impermeable  internal 
boundaries.  In  the  context  of  the  latter  topic  it  is  natural  to  include  the 
discussion  of  the  behaviour  of  the  bag  or  cloth  material  used  to  package  the 
propellant . 

2.3.1  Jump  Conditions  at  a  Discontinuity  in  Two  Phase  Flow 

1  0 

We  have  previously  derived  these  results  for  one  dimensional  flow  .  The 
extension  to  multidimensional  flow  is  straightforward^^  and  requires  no  new 
discussion.  We  designate  the  properties  on  each  side  of  the  discontinuity  by 
the  subscripts  1  and  2.  Moreover,  we  use  Uj^  and  Up  to  denote  the  normal  com¬ 
ponents  of  gas  and  solid  phase  velocity  and  Ur^-  and  Up^  to  denote  the  trans¬ 
verse  components.  The  jump  conditions  are: 


j  =  e.p.  (u  -  u  )  =  e.p.(u  -  u  ) 
s  2*^2  s 


2.3. 1.1 


j  =  (l-&,)p  (u  -  u  )  =  (l-£^)p  (u  -  u  ) 

P  P  S  ^  p  s 


2. 3. 1,2 


2  ^P  2 

p  -f  (l-e  )R  +  - (u  -  u  )  +  (l-£  ) — (u  -  u  )  2.3. 1.3 

^  K  1  g  ^  n  ^  rg  p  s 

o  1  o  n 
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^2^2  2  *^P  2 

p  +  (l-£  )R  +  - (u  -  u  )  +  (l-£  ) — (u  -  u  ) 

g  ^  n  2^g  ^  p 
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Pi  Wr  ,  P2  ^2-  ^s)^ 


2. 3. 1.4 
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2. 3*1. 5 


j  (u^  -  )  =  0 

^  1  2 


j 


P 


2.3. 1.6 


Here  we  have  used  Ug  to  denote  the  normal  velocity  of  a  point  on  the  surface 
of  discontinuity.  It  is  important  to  note  that  separate  momentum  jumps  are 
not  provided  for  the  individual  species.  The  momentum  jump  is  expressed  for 
the  mixture  alone.  In  general  an  exchange  of  momentum  between  the  phases  is 
not  forbidden  at  the  discontinuity. 

We  will  confine  our  interest  to  the  case  jp  -  0  since  we  do  not  consider 
boundaries  which  are  permeable  to  the  solid  phase.  We  may  then  distinguish 
two  cases  of  interest. 

(a)  Impermeable  Boundary,  j  =  0 

Evidently  the  boundary  conditions  are  the  intuitively  obvious  set.  We 

have : 


u  -  u  =  u  2.3. 1.7 

n  n  s 
1  2 


2. 3. 1.8 


+  (l-e2)R2 


2.3. 1.9 


and  the  tangential  components  Ur^^,  ,  Um  ,  u^  ,  u  are  unrestricted  corres- 

.  ^2  PTi  PT2 

ponding  to  slip  conditions.  In  the  event  that  the  flow  exists  on  only  one 
side  of  the  boundary  2. 3. 1.7,  2.3. 1.8  and  2.3. 1.9  apply  to  the  non-trivial 
side.  We  also  note  that  2.3. 1.7  is  the  familiar  condition  for  the  inviscid 
single  phase  flow. 

(b)  Gas  Permeable  Boundary^  j  /  0 

We  have  the  conditions: 
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2.3. 1. 10 
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2.3.1.13 


Also,  according  to  2. 3. 1.6,  an  arbitrary  jump  in  the  tangential  velocity  of 
the  solid  phase  is  admitted.  We  observe  that  as  ^2  ^  these  conditions 
reduce  to  the  familiar  Rankine-Hugoniot  conditions  for  a  single  phase  fluid. 
In  fact,  2,3.1.12,  the  condition  of  adiabatic  flow,  is  unchanged  by  the  pre¬ 
sence  of  the  solid  phase.  The  condition  of  continuity  of  mass,  2.3.1.10, 
appears  as  the  obvious  generalization.  The  major  novel  features  are  associ¬ 
ated  with  the  momentum  jump,  2.3,1.11.  We  see  that  progress  cannot  be  made 
without  either  specifying  the  manner  in  which  the  granular  stress  behaves  at 
the  discontinuity  or  providing  some  functionally  equivalent  datum. 

Since  we  have  assumed  R  to  embed  the  stresses  due  to  direct  contacts 
among  the  grains  it  is  natural  to  postulate: 


(1  (1  e2^®-2 


2.3.1.14 


whereupon  the  gas  phase  momentum  jump  becomes: 

£  1 P 1  9  ^2  p2 


p  _j - - -  (u  -  u  )  =  p  + 

'  So  ^  ^  ^  Sc 
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U  ) 

n  s 
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2.3.1.15 


However,  the  postulational  nature  of  2.3.1.14  should  be  emphasized.  The 
tensor  R  arises  formally  as  the  difference  between  the  average  stress  ten¬ 
ors  for  the  two  phases^  .  The  argument  that  it  embeds  only  the  influence 
of  granular  contact  is  plausible  only  at  low  relative  Mach  numbers.  At  high 
Mach  numbers  the  situation  is  complicated  by  the  influence  of  the  pressure 
distribution  around  each  particle  associated  with  drag. 


On  the  other  hand,  2.3.1.15  agrees  with  the  familiar  form  of  the  momen¬ 
tum  jump  used  to  analyze  flow  losses  at  a  sudden  increase  in  the  cross  sec¬ 
tion  of  the  f low  ’  of  an  incompressible  fluid. 
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2, 3, 2  External  Boundaries 


We  do  not  consider  the  possibility  of  flow  of  the  solid  phase  through  an 
external  boundary.  The  gas  permeable  boundary  is  of  interest  only  for  prob¬ 
lems  involving  poor  obturation  or  detailed  simulation  of  the  ignition  train. 
The  mass  transfer  relations  for  such  cases  can  be  deduced  according  to  stan¬ 
dard  isentropic  flow  theory^^and  are  not  discussed  here.  We  confine  our 
attention,  for  the  time  being,  to  the  fully  impermeable  boundary. 

Accordingly,  we  have  to  consider  simply  the  slip  boundary  conditions  ap¬ 
plied  to  both  phases,  equations  2.3. 1.7,  2.3. 1.8  and  2.3. 1.9.  For  the  sur¬ 
faces  of  the  tube  and  the  breech  we  assume  u  =0,  corresponding  to  the  ne¬ 
glect  of  recoil  and  strain. 

Now  let  Spj^Qj  be  the  surface  of  the  projectile  which  intrudes  into  the 
combustion  chamber  and  let  S  be  the  local  normal  to  the  surface,  positive 
pointing  outwards  from  the  combustion  chamber.  Let  n  =  n^)  describe  the 

2  -  and  r  -  components  of  S.  Let  the  axial  speed  of  the  projectile  be  Up^^^. 
Then  the  boundary  conditions  at  the  projectile  surface  are: 


u  • 


^z^PROJ 


2. 3. 2.1 


Of  course,  the  condition  on  the  solid  phase  applies  only  in  a  region  of  two 
phase  flow. 

The  projectile  motion  is  assumed  to  be  governed  by: 


T  r 

M  — — —  =  I  (P  +  (l-£)R)n  da  -  F  2.3.2 

eff  dt  I  ^  z  res 

PROJ 

where  the  effective  mass  is  related  to  the  actual  mass,  M,  and  polar  moment 
of  inertia,  I,  according  to: 


"eff  ■  "  +  5^ 


2. 3. 2. 3 


and  D  and  0  are  respectively  the  diameter  of  the  tube  and  the  angle  of  the 
rifling.  The  quantity  F^^^  is  the  bore  resistance  and  will  be  assumed  to 
be  available  as  a  semi -empirical  correlation. 
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2.3.3  Gas  Permeable  Internal  Boundaries 


We  take  equations  2.3.1.10,  2.3.1.12,  2.3.1.13,  2.3.1.14  and  2.3,1.15  to 
apply  at  all  the  internal  boundaries  defined  by  a  jump  in  porosity.  When  the 
boundary  separates  two  regions  of  two  dimensional  flow,  the  interpretation  of 
these  conditions  is  straightforward  since  they  refer  to  the  boundary  values 
on  each  side  of  the  discontinuity.  We  note,  of  course,  that  these  conditions 
reduce  correctly  when  one  or  both  sides  of  the  boundary  contains  a  single 
phase  flow. 


When  the  boundary  is  such  that  the  flow  on  one  or  both  sides  is  quasi- 
one-dimensional  or  lumped  parameter,  care  must  be  taken  with  the  interpreta¬ 
tion  of  terms.  For  example  let  side  1  correspond  to  two  dimensional  two  phase 


flow  so  that  values  of  e, 


and 


PTi 


are  available.  If 


1 

the  side  2  represents  a  quasi-one-dimensional  model  of  the  tangential  flow,  we 
identify  £2^  P2  with  the  state  variables  for  side  2.  However,  u^^ 
modeled  except  implicitly  through  the  mass  transfer  terms  and  m^.  The  values 
of  P- ,  e^,  uq^^  correspond  with  the  state  of  the  quasi-one-dimensional  flow  only 
when  that  region  is  acting  as  a  donor. 


The  assumption  that  2.3.1.14  and  2.3.1.15  are  valid  restricts  our  analy¬ 
sis  to  relatively  low  subsonic  transfers  through  the  internal  boundaries.  In 
our  opinion,  further  study  is  required  to  permit  a  meaningful  analysis  of  the 
transonic  and  supersonic  flows.  Certainly,  one  can  proceed  pragmatically, 
using  quasi-equilibrium  nozzle  flow  analyses  as  a  basis  for  determining  solu¬ 
tions.  However,  the  validity  of  such  an  approach  is  not  clear  at  present. 


Finally,  we  comment  that  from  the  application  of  these  boundary  condi¬ 
tions  we,  in  effect,  establish  the  constitutive  laws  for  the  mass  transfer 

terms  m.  and  m  as  used  in  section  2.1. 

1  o 

2.3.4  Impermeable  Internal  Boundaries 


This  topic  is  of  considerable  importance  since  it  is  here  that  we  embed 
the  influence  of  the  bag  and  additive  materials  on  mass  transfer  and  hence,  on 
flamespreading.  If  the  bag  is  impermeable  we  have  as  boundary  conditions, 

2.3. 1.7,  2.3. 1.8  and  2.3. 1.9.  As  in  the  previous  section,  the  interpretation 
of  these  terms  is  straightforward  provided  that  both  sides  of  the  boundary 
correspond  to  a  two  dimensional  flow. 

When  one  or  both  sides  is  a  quasi-one-dimensional  two  phase  flow,  the 
boundary  condition  2.3. 1.9  serves,  in  effect,  to  define  the  evolution  of  the 
cross  sectional  area.  As  we  have  shown  previously^^ ,  the  constitutive  law 
for  the  granular  stress  may  be  combined  with  the  continuity  equation  for  the 
solid  phase  to  produce  the  desired  result.  We  also  showed^  that  when  the 
bag  became  ruptured,  the  mass  transfer  deduced  according  to  the  previous 
section  could  be  used  to  estimate  the  subsequent  lateral  dilation  of  the 
quasi-one-dimensional  flow. 

If  the  data  base  warrants,  it  may  be  of  value  to  treat  the  internal 
boundary  b}^  explicitly  recognizing  the  flow  resistance  and  thermal  loss  as¬ 
sociated  with  the  bag  material.  Equations  2.3.1.11  and  2.3.1.12  may  be 
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modified  to  reflect  a  finite  loss  of  momentum  and  energy  experienced  by  gas 
passing  through  the  bag.  Indeed,  the  permeability  may  be  made  time  depend¬ 
ent  and  linked  to  flow  parameters  which  characterize  the  mechanical  and 
thermal  attack  due  to  the  penetrating  gas.  Moreover,  the  momentum  jump  may 
also  reflect  an  inertial  loss  associated  with  acceleration  of  the  bag  mater¬ 
ial.  This  level  of  modeling  may  well  be  necessary  in  order  to  study  charge 
configurations  composed  of  several  bags  or  even  for  single  bag  charges  in 
which  a  sophisticated  use  is  made  of  the  bag  permeability  to  control  flame¬ 
spreading. 
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2.4  Characteristic  Analysis  of  Balance  Equations 


Our  interest  in  this  topic  is  principally  motivated  by  the  numerical  ram¬ 
ifications  of  the  theory  of  characteristic  surfaces.  In  the  case  of  one  di¬ 
mensional  unsteady  flow,  the  existence  of  real  characteristic  directions  en¬ 
ables  one  to  replace  the  system  of  partial  differential  equations  by  an  equiv¬ 
alent  system  of  ordinary  differential  equations  in  which  the  derivatives  are 
taken  along  the  characteristic  lines.  When  we  proceed  to  a  larger  number  of 
independent  variables  an  analogous  result  holds  for  hyperbolic  systems  of  equa¬ 
tions.  Given  n  independent  variables,  a  hyperbolic  system  is  one  that  admits 
the  existence  of  a  hypersurface  of  dimension  n-1  such  that  only  derivatives 
interior  to  the  surface  appear  in  the  equations. 

We  proceed  as  follows.  In  section  2.4.1,  we  discuss  the  theory  in  gen¬ 
eral  for  a  quasi  linear  system  of  partial  differential  equations  which  depend 
on  three  independent  coordinates.  In  section  2.4.2,  we  illustrate  the  theory 
by  reference  to  the  equations  for  unsteady  homentropic  flow  with  azimuthal 
symmetry . 

In  section  2.4.3  and  2.4.4  we  deduce  the  characteristic  forms  for  two  di¬ 
mensional  single  phase  and  two  phase  flows  respectively. 

2.4.1  Formulation  of  Characteristic  Analysis 


Consider  a  system  of  partial  differential  equations 


+  B 


3z 


+  C 


9r 


2 . 4 .  1 . 1 


where  ip  and  D  are  n-dimensional  column  vectors  and  A,  B,  C  are  n  x  n  square  ma¬ 
trices,  The  concept  of  a  characteristic  surface  follows  naturally  from  the 
consideration  of  an  initial  value  problem  posed  for  a  surface: 


cj)(t,  z,  r)  =  0 


or,  in  general,  the  family  of  surfaces  generated  by  the  parameter  such 
that : 


ct)(t,  z,  r)  =  cj)^ 


2.4. 1.2 


Let  a,  3  coordinates  internal  to  the  surface  cj)  =  cj)  then  cp  itself  serves 
as  a  normal  coordinate.  On  cj)  =  (j)^  we  assume  that  we  are  given  values  of  ip 
and  hence,  values  of  ip^  and  ipo.  The  surface  (p  =  (p  is  said  to  be  free  if 
equation  2.4. 1.1  permits  the  determination  of  the  normal  derivative  and 
characteristic  if  it  does  not.  If  (J)  =  4)^  is  characteristic  it  follows  that 
2.4. 1.1  may  be  expressed  in  terms  of  derivatives  with  respect  to  a  and  3 
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alone,  that  is  to  say,  derivatives  internal  to  the  characteristic  surface. 


By  means  of  the  chain  rule  for  differentiation,  2.4. 1.1  may  be  trans¬ 
formed  into  derivatives  with  respect  to  (j),  a  and 


[A({)  +  Bcj)  +  C(})  =  D  -  [Aa  +  Ba  +  Ca  lil; 


[AB^+  B6^+  2.4. 1.3 


Accordingly,  the  question  of  whether  cj)  is  free  or  characteristic  is  settled 
by  the  rank  of  the  matrix: 


A  =  A(|)^  +  Bcj)^  +  Ccf)^ 


2.4. 1.4 


If  Rank  (A)  =  n,  the  system  2.4. 1.3  always  has  a  unique  solution  and  the 
surface  (J)  =  is  free.  However,  if  the  value  of  Rank  (A)  <n  then^2.4.1.3 
does  not  possess  a  solution  for  arbitrary  initial  data  on  the  surface  ([)  =  ({)q. 
In  the  latter  case  the  partial  differential  equation  represents  a  constraint 
on  the  data  as  expressed  by  the  condition  of  solvability  of  2.4. 1.3.  Thus  if 
we  let  A^  be  the  augmented  matrix  formed  by  appending  to  A  the  column  vector 
corresponding  to  the  right  hand  side  of  2.4. 1.3. 


A'*'  =  [A;  D  -  [Aa  +  Ba  +  Ca  -  [A3^+  B3  +  C3 


Then  the  condition  of  solvability  is 


57  . 


Rank  (A^)  =  Rank  (A) 


2.4. 1.5 


The  condition  Rank  (A) <  n  will  lead  to  a  partial  differential  equation 
for  (j)(t,z,r)  in  the  form: 


2.4. 1.6 


it  being  assumed  that  A,B,C  are  functions  only  of  t,z,r  and  if;.  Since  every 
element  of  A  is  a  homogeneous  linear  combination  of  (p^  ^  (p^  and  c})^  it  follows 
that  F  is  homogeneous  of  order  k  ^  1  in  these  quantities  so  that: 


F(t,z,r,i|j,X({)  ,A(])  .Atj)  ) 

L  Z  JU 


A'^F(t,z,r,(p,4)  ,(!)  ,([)  ) 
L  z  r 


^^■Hadley,  G.  "LindoA  klge-b/ui*'  AddUon-W^itzy  1961 
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Accordingly  it  follows  that; 


=  0 


2.4. 1.8 


Because  of  the  degree  of  freedom  induced  by  the  homogeneity  of  F  it  is 
convenient  in  many  cases  to  append  an  additional  condition  corresponding  to 
the  normalization  of  the  vector  (<Pt  9(t>z9‘Pr)  •  practice  the  most  convenient 

choice  is  to  set  (p^  =  -  !•  This  corresponds  to  having  (j)  in  the  form: 


Cj)  =  t(z,r)  -  t 


2.4. 1.9 


so  that  =  -  1,  (p^  =  dt/dz  and  =  dt/dr. 
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It  is  useful  to  interpret  these  results  geometrically  ,  We  may  think 
of  (})  =  (j)^  as  defining  a  surface  with  normal  vector  proportional  to  • 

Then  2.4. 1.6,  the  partial  differential  equation  for  the  characteristic  sur¬ 
face,  imposes  an  algebraic  constraint  on  the  components  of  the  normal  at  each 
point  in  the  (t,z,r)  space.  At  each  point  2.4. 1.6  defines  a  family  of  planes 
such  that  the  characteristic  surface  must  be  tangent  to  one  of  them.  If  F  is 
not  linear,  the  envelope  of  this  family  of  planes  is  a  cone,  the  Monge  cone, 
whose  generators  are  called  bicharacteristics. 

According  to  2,4. 1.7  the  family  of  allowable  normal  vectors  at  a  given 
point  lies  on  the  surface  of  a  cone  whose  apex  is  the  point  in  question.  Thus 
2.4. 1.8  asserts  that  the  vector  (F^j,^,  F^  ,  F^  )  is  perpendicular  not  only  to 
the  surface  of  normal  vectors,  but  also  to  the  vectors  themselves  since  they 
lie  along  the  cone.  In  fact,  the  vec tor  ( F(j)^ ,  Fcp^)  defines  a  bicharacter¬ 

istic  direction.  From  2.4. 1.8  it  is  evident  that  it  must  lie  in  a  tangent 
plane  of  a  charac-teristic  surface.  However,  the  bicharacteristic  may  be  thought 
of  as  the  limit  of  the  line  of  intersection  of  neighboring  tangent  planes.  Thus 
if  we  write  b  as  a  bicharacteristic,  n  as  the  normal  to  a  tangent  plane  and  A  as 
a  parameter  which  labels  the  planes  at  a  given  point,  it  fallows  that  b  =  n  x 

(n  4-  --  dA).  Therefore  b  =  n  x  dA .  Since  both  n  and  ~  lie  on  a  tangent 
dA  nA  aA 

plane  of  the  cone  defined  by  2.4.  1.6  we  see  that  b  is  parallel  to  (1^^  ,  F^  , 

Fa  ).  The  bicharacteristic  ray  may  be  written  as:  t  z 


dt  dz  dr 


2.4.1.10 


^^•CouAant,  R.  and  H-lCbeM,  V.  ''MeXhod-6  Mathmat-icai  Phiplci"  Inte^-iclence. 
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This  result  may  be  used  to  eliminate  (p^  and  (j)^  from  2.4. 1.6  and  to 
describe  the  characteristic  surface  by  reference  to  the  bicharacteristics. 

2.4.2  Illustration;  Homentropic  Flow 

In  order  to  understand  the  formulism  we  consider  briefly  the  balance 
equations  for  homentropic,  unsteady,  single  phase,  compressible  flow  in  two 
dimensions  with  cylindrical  symmetry. 


9t 


+ 


A 

9z 


pu  + 


_9_ 

9r 


pv 


r 


0 


2.4. 2.1 


9u 

9t 


+  u 


V 


9u 

9r 


9P  - 


p  9z 


=  0 


2. 4. 2. 2 


u 


9v 

9z 


+  V 


dr 


=  0 


2.4. 2.3 


Here  we  have  used  c^  =  (9p/9p) 


s 


as  the  square  of  the  speed  of  sound. 


We  have; 


V 

— 

-  p  - 

p 

r 

u 

D  = 

0 

V 

0 

A  = 


lo 

o 

LI 

u  p  0 

0  1  0 

;  B  = 

c^/p  u  0 

o 

o 

_j 

0  0  u 

V  0  p 

C  = 

0  V  0 

c^/p  0  V 
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Thus  we  must  consider  the  rank  of: 


+  U(})^  +  v4)^ 


-  -  -  ij, 

P  z 


p  'r 


4>t  +  0 


0  d)  +  ud)  +  vd) 

z 


A  necessary  condition  that  the  rank  of  A  be  less  than  3  is  det(A)  =  0  or: 


((})_^+  U(j)^+ 


C^((p  ^-+  (p  ^)}  =  0 


2.4. 2.4 


Evidently  the  characteristic  surface  consists  of  tv/o  sheets: 


=  (})^  +  u(j)^  +  v4)^  =  0 


2. 4. 2. 5 


and 


=  (‘J’t  ■*'  ~  °  2. 4. 2. 6 


The  bicharacteristics  corresponding  to  2.4. 2.5  are  just  the  particle  path¬ 
lines.  Those  for  the  quadratic  sheet  may  be  expressed  as: 


b  =  {2(ci)^  +  u(J)^  +  vcf)^),  2u((})^  +  u(t)^  +  v(f)^)  -  2c2(j)^, 

2v(d)  +  ucb  +  vd)  )  -  Zc^cp  } 
t  z  r  r 


Accordingly  we  can  express  (j)^,  (p^  and  (p)^  in  terms  of  the  components  of 
b  =  (dt,  dz,  dr)  as: 


(p)_^  =  {dt(c^“  u^-  v^-)  +  udz  +  vdr}/c2 
d)  =  (udt  -  dz)/c^ 


d)  =  (vdt  -  dr)/c^ 
r 
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whereupon  substitution  into  2.4. 2,6  shows  the  bicharacteristics  to  satisfy: 


-  U)^  +  (f  -  V)^  = 


2.4. 2. 7 


In  a  computational  application  we  expect  that  dt  is  known  as  the  time 
increment  through  which  the  solution  is  to  be  advanced.  Moreover,  we  select 
the  bicharacteristic  to  lie  in  a  given  plane.  Therefore  the  ratio  dz:dr  is 
known  and  2,4.2. 7  enables  us  to  determine  their  separate  values.  These,  when 
entered  into  the  condition  of  compatibility  lead  to  a  finite  difference  formula 
for  the  state  variables.  Since  we  have  no  application  for  this  simple  case  we 
do  not  bother  to  write  down  the  condition  of  compatibility.  It  follows  mechan¬ 
ically  from  a  consideration  of  the  augmented  matrix. 

2.4.3  Single  Phase  Two  Dimensional  Unsteady  Flow 

We  now  wish  to  deduce  complete  results  for  a  two  dimensional  single 
phase  flow.  These  results  will  differ  from  those  of  the  preceding  section 
in  three  essential  respects.  We  will  assume  that  an  initial  transformation 
of  variables  has  been  made  for  computational  purposes.  Secondly,  we  will  not 
take  the  entropy  to  be  the  same  for  all  fluid  particles.  Thirdly,  we  will  de¬ 
duce  the  conditions  of  compatibility. 


The  balance  equations  are  given  in  cylindrical  coordinates  in  section 
2.1.3.  As  usual,  we  recast  the  energy  equation  by  means  of  the  continuity 
equation  and  the  identity: 
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3P. 


g 


2.4.3. 1 


Thus  we  have: 
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2.4. 3.2 
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o  9z 


-  =  Kz 


2.4. 3. 3 
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2.4. 3.4 
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where  ip  is  understood  to  be  a  source  term  and  is  not  to  be  confused  with  the 
column  vector  of  state  variables. 

We  identify  these  equations  with  the  system: 

A-|^+b|^-+c|^=  D  2. A. 3. 6 

dt  02  dr 


by  setting: 


At  this  point  we  have  established  the  balance  equations  in  a  form  suit" 
able  for  the  application  of  the  methodology  described  in  section  2.4.1.  How¬ 
ever,  we  now  consider  a  transformation  of  coordinates  in  the  form: 

T  t 

=  c(t,z,r)  2.4. 3.8 

n  =  n(t,z,r) 

We  assume  that  this  transformation  is  one  to  one  and  has  continuous  partial 
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derivatives  and  that  9(C,ri)/9(z5r)/"  0  so  that  we  can  also  write: 


t  =  T 

2  =  z(T,C,ri) 

r  -  r(T,^,n) 


2.4. 3.9 


We  have  used  a  separate  notation  t  for  the  time  in  the  transformed  frame.  This 
facilitates  the  use  of  subscripts  to  denote  partial  derivatives.  Thus  we 

understand  cj)  =  ('^)  whereas  ^  =  (-^)  and  cj)  is  an  arbitrary  property. 

■t  z,r  C>n 

Bearing  this  in  mind  we  define: 


^  2.4.3.10 

V  =  r 
m  T 

Thus  u^  and  v^  are  the  velocity  components  in  a  cylindrical  coordinate  frame 
of  a  point  moving  so  that  it  is  stationary  with  respect  to  the  transformed 
frame.  Evidently,  if  we  impose  the  requirement  Uj^  =  u  and  v^  =  v  we  will 
have  selected  the  transformed  frame  to  coincide  with  a  Lagrangian  description 
of  the  fluid  whereas  the  choice  u^  =  v^  =  0  implies  the  retention  of  an 
Eulerian  description,  possibly  in  a  different  coordinate  frame  established  by 
a  stationary  transformation. 

It  follows  that  the  balance  equations  2. 4. 3. 6  subject  to  2.4. 3.8  and 
2.3.4.10  become: 

aM  +  [(B-Au^)C/  (C-Av^)?^]||  +  [(B-Au^)n^+  (C-Av^)n^]  |^=  D 

2.4.3.11 


Thus  we  now  consider  the  characteristic  surfaces  for  the  system: 


where  we  identify  B’  and  C’  as: 
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and  where  we  have  introduced: 


W  =  (u-u^)52  + 


X  =  (u-u^)n^  +  (v-v^)n^ 


2.4.3. 13 


Then  the  characteristic  surfaces  c})(T,^,n)  such  that  the  rank  of  A 

is  less  than  four  where; 


A  -  Ad)  +  B’d)  +  C’d) 

'n 


2.4.3. 14 


Thus  we  have: 


C^cJj/g 


p(^  (p  +  r}  (p  ) 


g  (C^(p^  +  n  fp  ) 

o  z  <;  z  ri 

g  (?  (()  +  ri  A) 

o  r  ^  r  n 


and  where  we  have  introduced: 


=  d)  +  W(l)  +  xcf) 

’  T  ^  n 


2.4.3.15 


Now  in  our  applications  of  the  conditions  of  compatibility  we  shall 
always  require  that  either  <()^  =  0  or  that  (f)^  =  0  so  that  the  normal  lies 
either  in  the  T-r]  plane  or  the  T-C  plane.  We  assume  therefore  that  (})^  =  0. 
The  corresponding  results  for  the  case  (j)^  =  0  will  follow  from  considera¬ 
tions  of  symmetry.  With  this  assumption  A  reduces  to: 


PC  4> 
r  ? 


o  z  ^ 


2.4.3.16 
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Now  we  consider  two  possibilities 
(i)  Let  (^  -  0  Then  A  reduces  to: 


From  this  it  is  apparent  that  the  streamline 


^  =  (j)^  +  wc})^  =  0 

is  a  characteristic  direction.  In  order  to  establish  the  condition  of  com¬ 
patibility  we  now  introduce  a  and  3  as  coordinates  internal  to  the  charac¬ 
teristic  surface.  Then  we  may  write  the  augmented  matrix  as: 


{ap  +  a  +  n  ot  lu  +  pr^  a  +  ri  a  Iv  } 
'  z  C  z  n  a  r;  'r  a 


{3P3  +  +  n^3Ju^  +  P[C/^  +  \3Jv^} 


5  -  {pav  +  g[i;a  +rict]p} 

3  a  ^o'-  r  r  it'  a 


-  {p8vg  + 


\  -i[P„  -•P"P„/8j-  BfPg  -  c2p^/gj 


where  a  and  3  are  defined  by  analogy  with  (j. 
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Accordingly,  if  we  write  i=l,...,  4  to  denote  the  members  of  the 
fifth  column  of  A'*’,  the  conditions  of  solvability  yield  the  following  con¬ 
ditions  of  compatibility. 


5  -  0 

4 


2.4.3.17 


2.4.3. 18 


A  convenient  choice  for  a  ,  3  is: 


a  =  C 


3  =  n 


2.4.3.19 


and  we  will  adhere  to  this  convention.  The  use  of  a  separate  nomenclature 
for  the  coordinates  internal  to  the  characteristic  surface  is  again  motivated 
by  the  desire  to  maintain  clarity  in  respect  to  the  representation  of  the  par¬ 
tial  derivatives. 

Evidently: 


a  =  a 
T  n 


=  0 


a  =  3  =  1 

4  n 


It  follows  that  wv^  =  v^  +  wijj^  whereupon  2.4.3,17  may  be  identified  as  a 
linear  combination  of  the  two  momentum  equations. 

On  the  other  hand,  since  the  bicharacteristic  ray  satisfies  da  =  wdT  we 
can  write  2.4.3.18  as: 


[p - P  ]da  =  [C  -  (p  -  —  Pfl)]dT  2.4.3.20 

'■‘^a  g  01^  ‘-^4  g  ■' 

o  o 

This  is  now  recognized  as  the  familiar  one  dimensional  result  with  the 
3  -  derivatives  (i.e.  q  -derivatives)  taken  to  the  right  hand  side  and  treated 
formally  as  non-homogeneous  terms. 

(ii)  Now  let  (()  ^  0.  Then  perform  successively  the  following  row  and  column 
operations  to  A  as  given  by ^equation  2.4,3.16.  Add  times  column  one  to 

column  four;  subtract  times  row  two  from  row  one;  subtract  g^^2(})^/p(f 

times  column  two  from  column  four;  subtract  g^^^cp^/p^  times  column  three  from 
column  four;  subtract  -  times  row  three  from  row  one;  add  go/c^  times  row 

four  to  row  one.  Then  A  is  equivalent  to  A  where: 
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0  0  [^2_  c2(r  2(j^  2  +  ;  2^  2)] 


(f)C 


0 


pc()  0 

0 


0 

-c^4)/g^  0  0 

Accordingly,  (f)  can  only  be  characteristic  if: 


cj>^  = 


That  is: 


2 

0 

0 


(})  =  -[w  +  c/c  2+  ^  2]^ 

f-  '7.  r>  4 


z  r  “  ■  ^ 

Thus  the  bicharacteristics  satisfy  the  familiar  one  dimensional  form: 


dx  = 


w  +  c(5 


The  corresponding  condition  of  compatibility  is  easily  seen  to  be: 

1  9  2  9  3^4 

Then  choosing  a  and  3  before,  equation  2.4.3.19,  and  observing  9^/9 
where  c^^  =  condition  of  compatibility  may  be  expressec 
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As  before,  the  3  ”  derivatives  are,  in  effect,  n  derivatives  which  is  to 
say  derivatives  along  a  coordinate  curve  which  we  may  align  with  a  compu¬ 
tational  boundary.  Again,  2.4.3.24  is  analogous  with  the  one  dimensional 
result  with  the  cross  derivatives  (3  -  derivatives)  treated  formally  as  non 
homogeneous  terms. 

2.4.4.  Two  Phase  Two  Dimensional  Unsteady  Flow 

Finally,  we  wish  to  establish  complete  results  for  a  two  dimensional 
two  phase  flow,  including  the  conditions  of  compatibility.  As  in  the  pre¬ 
ceding  section  we  shall  see  that  the  algebraic  problem  is  complicated  only 
to  the  extent  that  the  non-homogeneous  terms  now  include  derivatives  along 
one  of  the  internal  coordinate  directions.  This  direction  will  always  con¬ 
form  with  a  computational  coordinate  line  as  a  consequence  of  the  transfor¬ 
mation  2.4. 3.8  which  we  will  also  apply  to  the  two  phase  flow  equations. 

In  addition,  we  will  introduce  a  parameter  A  into  the  solid  phase  mo¬ 
mentum  equations  such  that  when  A  =  1  the  pressure  gradient  is  treated, 
correctly,  as  a  differential  term  and  when  A  =  0,  the  pressure  gradient  is 
treated  as  a  non-homogeneous  term.  It  is  emphasized  that  the  results 
corresponding  to  a  non-unit  value  of  A  are  to  be  thought  of  as  merely  for¬ 
mal;  they  do  not  represent  the  hyperbolic  structure  of  the  balance  equa¬ 
tions.  Our  purpose  in  introducing  them  is  to  provide  certain  approximate 
results  for  use  in  determining  numerical  solutions. 

The  system  to  be  studied  includes,  of  course,  the  balance  equations 
2.1.1,!  through  2. 1.1. 7.  In  addition,  we  must  consider  the  constitutive  law 
2. 2. 2. 2  which  is  in  differential  form.  The  constitutive  laws  2. 2, 6. 2  and 
2.2.7. 1  need  not  be  considered.  Although  they  are  expressed  as  partial 
differential  equations,  they  are  not  coupled,  at  a  differential  level,  to 
the  other  governing  equations. 

Eliminating  the  internal  energy  in  favor  of  the  pressure  and  density  in 
2. 1.1.5  and  adding  the  appropriate  component  of  (A-l)  (l“e)gQVp  to  each  side 
of  2. 1.1. 7  leads  to  the  system,  in  transformed  coordinates: 


A  +  B'  +  C  =  D 

dT  9^  9ri 


2.4.4. 1 
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ex  epn^  ePHj, 


0  px 


0  epx 


0  eg  n  0 

o  z 


epx  eg  n  0 
or 


2.4.4. 6 


X  0 


0  Xp  -(l-e)n2  ° 


0  0  (l-e)PpXp  0  g^n^ 


0  A(l-e)g^n^  0 


0  (l-e).p  X  g  n 

p  p  o  r 


0  p  ^  0 


As  in  the  previous  section  we  have  introduced  the  kinematic  coefficients 

w,  X,  w  and  x  such  that: 

P  P 

w  =  (u~u  +  (V“V 

m  z  m  r 


X  =  (u-u  )n  +  (v-v  )n 
m  z  m  r 


W  =(U"U)^  +(v-v)^ 

p  p  m  z  p  m  r 


X  =  (u  -  u  )n  +  (v  -  V  )n 

p  p  m  z  p  m  r 


2.4.4. 7 


The  characteristics  are  determined  from  the  rank  of  A  given  by: 
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£(}) 


0  ep(|) 


+  n  (jj  ) 

r  n 


•"  Iz'l'n^ 


+ 


0  g^i(l-e)(c^<f^  0 

+  T1  ({)  ) 
z^n 


■(1-£)(C. 


(l-e)(.) 


0 


•(l-e)(C^({>^^  0 

+  Vn> 


0  2  ( r  d) 

+  n  j)  ) 
z  n 


0  0 

+  n  4>  ) 


+ 


0  p 

P 


2.4.4, 8 


and  where 


d)  =  d)  +  wd)  +  xd) 
^4  n 


h  =  d)  +wd)  +  xd) 
P  P  4  P^h 


We  do  not  pursue  the  general  case  any  further.  Our  interest  will  be  confined  to 
cases  in  which  (j)^  or  (|)^  =  0.  As  in  the  previous  section,  we  proceed  under  the  as¬ 
sumption  (j)^  =  0.  Then  the  characteristic  conditions  are  found  to  be: 


2  2  22  Ip  Pp2  2 

[(  —  +  w  )  -  c..  ]  [(—  +  Wp)  -  a.  ]  =  A(— )  — «  [--  +  w] 

4  4  P  4 


2. 4. 4. 9 
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2.4.4.10 


(})  =  0 
P 


2.4.4.11 


which  are  formally  identical  with  the  one  dimensional  results  .  We  have  used: 


2  +  ^2) 


2  =  .2/.  2 


z  r 


It  is  also  clear  why  we  have  introduced  the  coefficient  X,  When  X  =  0,  equa¬ 
tion  2,4.4. 9  factors  and  always  has  four  real  roots. 

With  the  same  choice  of  coordinates  internal  to  the  characteristic  sur¬ 
face  as  we  used  in  the  previous  section  we  may  tabulate  the  members  of  the  last 
column  of  the  aug.mented  matrix,  i=lj  8  as: 


5  ’  =  5  -  [ewp  +  epr  u  +  epr  v  +  pwe  ] 

1  I  *'  a  z  a  r  a  a-' 

-  [exp^  +  epn^u^  +  +  P^p] 

5  '  =  C  -  [epwu  +  eg  i;  p  ]  -  [epxu  +  eg  n  p.] 

2  2^  o  z  a  ‘•  3  o  z  3 

K'  =  i  -  [epwv  +  eg  ?  P  ]  -  [epxv  +  eg  p  p^] 

3  3  a  or  a  3  o  r  3 

c2  c2 

E  ^  =  r  -  w[p  -  —  p  1  -  xfp  -  —  p  1 

5  '  =  ?  -[we  -  (1-e)^  u  -  (l-elc  v  ] 
55pa  zpa  r  pa 

-  [x  e_  -  (l-e)ri  u  -  (l-e)Ti  v  ] 

p  6  z  p3  r  pe-" 

5  '  =  5  -  [X(l-e)g  ?  P  +  (l-£)p  w  u  +  g  ^  a  ] 

0  6  o  z  a  P  P  pet  o  z  a 

-  [A(l-e)g  n  Po  +  (l-e)p  X  u  +  g  ^  0-] 

o  z  3  "^p  p  p3  o  z  3 

5  5  -  [A(l-e)g  C  P  +  (l-e)p  w  V  +  g  c  a  ] 

7  7  ora  p  p  pa  o  r  a^ 

-  [A(l-e)g  n  P„  +  (l-£)p  X  V  +  g  p  a„] 

*■  o  r  3  P  P  P3  o  r  3 


C  '  =  5  -  w  [a  +  —  e  ]  “  X  [o„  +  —  £„' 

8  8  P  “  8  a  p  3  g  3' 

o  o 
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In  terms  of  these  quantities  we  may  express  the  conditions  of  compatibility 
corresponding  to  2.4.4.10  and  2.4.4.11  in  the  simple  forms: 


C  '  =  0 

4 


2.4.4.12 


5  '  =  0 
8 


2.4.4.13 


The  condition  of  compatibility  corresponding  to  the  acoustic  characteris¬ 
tics  of  equation  2.4.4.10  has  the  form: 


8  C. 


5  P 


P^P 

pf'- 


(5  '  -  5  0 

7  8 

+  (C  -  - 

O  ^5  Y  y. 

o  z  C 

C  ’)] 

Pp$P  S 

-e- 

■o 

8 

^  eg 

[?  ’  +  -T 

-  ^  I - -  — S-  ^  t 

C 

• 

C  ’  ] 

(t)p  1  c 

4  (j)  2 

3 

2.4.4.14 


In  order  to  deduce  results  in  a  form  suitable  for  our  computational 
applications,  we  define: 


5  ”  [exp  +  epp  u  4*  epp  v  +  pxe  1 

1  n  zp  rp  p"' 


5  -  [epxu  -I-  eg  n.,P„] 

2  P  o  z  p 


C  -  [epxv  +  eg  p  p 
3  n  o  r  p- 


5  -  x[p - p  ] 

4  n  8^  n-* 


5  “  (l-e)p  u  -  (l-e)p  V  ] 

5  P  ^  z  pp  .  r  pp 


5  -  [A(l“£)g  p  p  -f-  (l-c)p  XU  +  g  p  a  ’ 

6  o  z  p  P  P  pn  o  z  P' 


C  ~  [A(l-e)g  p  p  +  (l-c)p  XV  +  g  n  cr  1 
7  o  r  p  ^p  p  pp  r  p-^ 


p 

5  -  X  [a  +  —  e  ] 

8  P  n  8^  p-* 
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Next,  we  use  2.4.1.10  to  recast  2.4.4. 9  as: 


2  2 
[(^  -  w)  - 


- 


1 


•'  *  p  £  dx 

p 


2.4.4.15 


Evidently,  when  X  =  0  this  yields  the  roots: 


4^  =  w  +  c .  and  4^  =  w  +  a. 
dx  —  *  dx  D  —  * 


We  introduce 


di; 

y  =  "Ttr  -  w 
■’  dx 


2.4.4.16 


4$. 

y  =  j  “  w 

•^D  dx  r 


2.4.4.17 


and,  by  reference  to  2.4.4.14  we  define: 


eg  C  ? 

^ ^  ^  <^2  ^ .  -^7  ^  _  TT  ^  . 


c2 


y  2 


2.4.4.18 


*  -  ®o^r 


-  (5  + 

p  y  7  y 

p  p  '  ] 


5  ) 


Q  Z 


r  cr  r 

— ^  (?  + 

p  y  f;  y 

p  p  ^  ] 


5  ) 
8 


2.4.4.19 


It  then  follows,  when  X  -  0  and  cj)^  =  0  and  ~  ^  ^^9  corresponding  to  the  solid 
phase  acoustic  characteristic,  that  2,4.4.14  yields: 


a 


p  y 

-P  P 


S00+ 

p 


w 


{■ 


)  (1  +  -^)  (C„u  ,  +  C  V  _  ) 


p  z  pg 


4  2  +  r  2 


.  pg  -  } 


2.4.4.20 
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Moreover,  when  we  have  simply  (J)^  =  0  then  2,4.4. 14  may  be  written  as. 


yp  = 


“  1  +  —  (1  +  ^)(y  a..2) 


y  'p 

^p 

pyy  (l-e)(l  +3r“)(?zV 

_ p__ _ 

(1  + 


+  ep(<;  u  +  ^  V  )  + 
z  a  r  a 


g  p  w 

“p~  y  +  3^)  +  c/) 


a 


2  -  a..2) 


+  py 


e 

a 


2.4.4.21 


We  note  that  whereas  2.4.4.20  requires  A  =  0,  2.4.4.21  does  not.  In  practice, 
however,  we  will  use  A  =  0.  Evidently  2.4.4.21  reduces  to  2.4.4.20  when 
yp  =  +  a^,^.  We  will  use  2.4.4.20  to  determine  values  of  a  on  the  computational 
boundaries.  Equation  2.4.4.21  will  be  used,  with  y  =  +  c,  corresponding 
to  the  gas  phase  acoustic  characteristics,  to  determine  values  of  p.  We  also 
note : 


a 

a 


a  /(f) 

T  dx 


2.4.4.22 


All  the  foregoing  results  were  obtained  subject  to  the  assumption  c})^  =  0. 
The  corresponding  results  for  ~  0  follow  from  considerations  of  symmetry. 
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3.0  TWO  DIMENSIONAL  CONVECTIVE  FLAME SPREADING  IN  A  CLOSED  CHAMBER 


From  a  consideration  of  the  systems  of  balance  equations  presented  in 
the  preceding  section  it  may  be  seen  that  a  logical  step  in  the  development 
of  a  fully  two  dimensional  interior  ballistic  code  is  the  establishment  of  a 
technique  for  the  numerical  solution  of  the  equations  of  two  phase  flow  in  a 
two  dimensional,  axis3mmietric  chamber.  Techniques  have  been  developed  pre¬ 
viously  for  one  dimensional  two  phase  flows  in  the  full  degree  of  generality 
required  by  our  ultimate  goal.  In  the  present  section  we  confine  ourselves 
to  the  problem  of  convective  flamespreading  in  a  closed  chamber,  without 
ullage.  Thus  we  are  concerned,  for  the  moment,  only  with  impermeable  sta¬ 
tionary  external  boundaries  and  not  at  all  with  internal  boundaries.  How¬ 
ever,  we  will  consider  the  complications  imposed  by  the  detailed  geometry  of 
typical  gun  chambers,  namely  the  hemispherical  breech  closure  plug  and  the 
projectile  intrusion  due  to  the  presence  of  a  boattail.  The  manner  of  ex¬ 
tension  of  the  present  techniques,  in  combination  with  those  previously  de¬ 
veloped,  will  be  discussed  in  chapter  4.0. 

The  method  of  solution  for  a  single  region  of  continuous  two  phase  flow 
bounded  by  stationary  impermeable  walls  is  contained  in  section  3.1.  Solu¬ 
tions  generated  with  this  method  are  presented  in  section  3.2. 

3. 1  Method  of  Integration 


The  present  section  has  five  subsections,  reflecting  five  aspects  of 
the  task  of  numerical  integration.  In  3.1.1  we  discuss  the  choice  of  co¬ 
ordinate  transformation  used  to  map  the  physical  domain  onto  a  unit  square. 
In  3.1.2  we  discuss  the  method  of  integration  for  interior  mesh  points  and 
in  3.1.3  we  discuss  the  treatment  of  points  which  lie  on  the  boundaries.  In 
section  3.1.4  we  comment  on  the  special  treatment  reserved  for  the  boundary 
points  which  lie  on  the  corners  of  the  computational  domain.  Finally,  in 
section  3.1.5  we  note  some  numerical  devices. 

3.1.1  Coordinate  Transformations 


The  physical  domain  has  stationary  boundaries  in  the  present  study. 
Accordingly,  a  single  transformation  is  sufficient  to  map  the  physical  do¬ 
main  once  and  for  all  onto  a  unit  square.  However,  in  subsequent  work  we 
will  be  confronted  with  the  more  general  problem  in  which  the  boundaries  are 
mobile  and  which  therefore  demands  a  time  dependent  mapping  technique. 

The  approach  taken  here  consists  of  two  steps.  Initially,  we  transform 
the  arbitrarily  configured  physical  domain  onto  a  square  by  means  of  an 
equipotential  map.  Subsequently,  we  effect  a  solid  phase  Lagrangian  trans¬ 
formation  so  that  the  mesh  moves  with  the  particles.  This  second  step,  which 
is  not  really  required  in  the  present  study,  is  expected  to  provide  the  de¬ 
sired  extension  to  the  case  in  which  the  boundaries  are  mobile.  It  has, 
moreover,  several  benefits  even  in  the  present  case.  The  troublesome  convec¬ 
tive  terms  in  the  solid  phase  balance  equations^ ^  vanish.  Accordingly,  the 
path  of  flaraespreading  is  more  sharply  resolved  due  to  the  elimination  of 
numerical  diffusion  in  the  solid  phase  thermal  equation,  2. 2. 6. 2. 
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Our  approach  to  the  first  part  of  the  coordinate  transformation  follows 
that  of  Thompson  et  al^^ .  The  computational  coordinates  are  embedded  into 
an  elliptic  equation  as  follows: 


It  should  be  noted  that  as  we  have  neglected  the  divergence  terms,  3. 1.1.1 
and  3. 1.1.2  do  not  express  Laplace’s  equation  in  cylindrical  coordinates. 

We  do  not  solve  3. 1.1.1  and  3. 1.1.2  directly.  Rather,  we  solve  the 
inverted  system: 


3. 1.1.3 

3. 1.1.4 


6  =  z  z  +  r  r 

4  n  c  n 

Y  =  z  ^  +  r  ^ 

Thompson  et  al  solved  3. 1.1.3  and  3. 1.1.4  subject  to  Dirichlet  data  on  all 
boundaries.  In  effect,  therefore,  only  the  distribution  of  mesh  points  on 
the  boundaries  is  required.  The  use  of  3. 1.1.3  and  3. 1.1.4  produces  a 
smoothly  varying  network  of  mesh  lines  on  the  interior.  In  order  to  solve 
3. 1.1.3  and  3. 1.1.4,  we  replace  all  derivatives  by  second  order  differences. 
The  resulting  finite  difference  equations  are  solved  using  the  method  of  suc¬ 
cessive  over-relaxation^^. 

In  the  present  work,  Dirichlet  data  are  generated  on  each  boundary  by 
specifying  the  element  as  a  string  of  straight  line  segments.  Mesh  points 
are  always  located  on  the  ends  of  each  segment.  The  interior  of  each  seg¬ 
ment  is  assigned  mesh  points  at  regularly  spaced  intervals.  The  number  of 
points  assigned  to  the  interior  of  each  segment  may  be  prespecified.  If 
the  total  number  of  points  allocated  to  the  boundary  element  is  not  exhausted 
by  the  foregoing  criteria,  the  surplus  is  allocated  so  as  to  minimize  the 
maximum  physical  mesh  interval  in  the  given  boundary  element. 

1972 
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In  addition,  we  have  encoded  the  capacity  to  apply  Neumann  data  on  any 
or  all  boundary  segments  to  express  the  requirement  that  the  mesh  lines  in¬ 
tercept  the  boundaries  at  right  angles.  Each  of  the  four  boundary  elements 
may  only  have  one  of  Dirichlet  or  Neumann  data  over  each  of  its  segments. 
Otherwise,  the  selection  is  arbitrary.  When  Neumann  data  are  required  only 
the  segment  end  points  are  predetermined.  Boundary  values  are  chosen  to  sat 
isfy  the  orthogonality  condition  as  expressed  by  a  second  order  one-sided 
difference  formula  after  each  sweep  of  the  interior  by  the  relaxation  scheme 

3.1.2  Integration  at  Interior  Mesh  Points 


The  integration  scheme  used  for  interior  mesh  points  is  basically  the 
MacCormack  algorithm^^  with  some  refinements  suggested  by  Moretti^^.  The 
basic  algorithm  of  MacCormack,  applied  to  the  system: 


If  +  B  +  C  =  D 

dt  9^  9ri 


3. 1.2.1 


may  be  expressed  as  a  two-level  scheme: 


^  +  [D.^  -4^  . 

1+1,3 


I  \ 

I ' .  . ) 
1,3 


C.'^. 

(ip.”-.  .  -  ip.^.  )]  At 

An  '^1,3+1  ^1,3 


3. 1.2.2 


-  'it’- 


c.  .  ^ 

1  5  /  I  ,  \  At 

—  »i.3  -n.j-i)!  T 


3. 1.2.3 


Thus  alternating  forward  and  backward  differences  are  used  for  the  repre¬ 
sentation  of  the  spacewise  derivatives. 

The  modification  suggested  by  Moretti  relates  to  the  discretization 
of  the  convective  derivatives.  These  are  always  represented  by  upstream 
differences  as  follows,  except  where  forbidden  by  proximity  to  a  boundary: 


Predictor 


9  _  +  1  r  . 

9c  “  Ac 


([).  .] 


3. 1.2.4 


Corrector 


9(1)  _  4-  1 
9C  "  AC 


[3(}) 


i+l,j 


^i+2,j 


3. 1.2.5 
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The  upper  or  lower  sign  is  used  according  as  the  pre-multiplying  velocity 
component  is  negative  or  positive  respectively.  It  should  be  noted  that 
3. 1.2.5  is  not  a  second  order  accurate  form.  It  only  yields  formal  second 
order  accuracy  in  combination  with  3. 1.2.4.  When  the  mesh  point  is  adjacent 
to  a  boundary  and  the  rule  expressed  by  3. 1.2.4,  3. 1.2.5  would  require  data 
outside  the  computational  domain,  we  revert  to  the  regular  MacCormack  pre¬ 
scription. 


We  also  note  that  the  matrices  B  and  C  involve  terms  like  These 

are  deduced  by  first  expressing  z^,  z^ ,  r^,  r^  by  means  of  centered  differ¬ 
ences.  Then  follow  from  an  elementary  theorem  of  partial 

differentiation  as; 


3. 1.2.6 


where  J  =  z^r  -  z  r^. 

?  n  n  ? 

It  was  not  found  necessary,  in  the  cases  considered  in  the  present 
study,  either  to  time-split  the  non-homogeneous  terms  or  to  introduce  im¬ 
plicitness  for  stability.  However,  we  have  incorporated  an  element  of  im¬ 
plicitness  into  the  determination  of  the  interphase  drag.  If  we  write 
fc.  =  (u-u^)  cj)  (u-Up.)  then  the  vectorial  prefactor  is  represented  implicitly 
in  the  term  u. 


The  integration  scheme  is  assumed  to  be  stable  when  subjected  to  a 
usual  C-F-L  domain  of  dependence  limitation.  If  C  is  the  fastest  local 
wavespeed  we  require: 


C 


^  < 
A?  - 


(r  -kr  )z  +  (z  -kz  )r 

n  c  4  n  c  ^ 


/(r  -kr  )  2  +  (z  -kz 

n  ^  n  C 


3. 1.2.7 


Where  k  =  +  according  as  practice,  we  further 

constrain  this  heuristic  limit  by  dividing  it  by  a  safety  factor  which  we 
have  taken  to  be  1.1. 


3.1.3  Integration  at  Boundary  Points 

31 

At  the  boundary  points  we  use  a  variation  by  Moretti  of  a  scheme 
apparently  due  to  Kentzer'^  .  The  method  is  based  on  the  characteristic 
forms  of  the  balance  equations.  However,  the  characteristic  forms  are  in¬ 
tegrated  explicitly  using  one  sided  differences  for  the  spacewise  deriva¬ 
tives.  We  describe  first  the  computational  sequence  and  subsequently  the 
procedure  of  discretization. 
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Along  a  C  =  constant  boundary  we  may  combine  the  components  of  the  gas 
phase  momentum  equation  to  yield  a  tangential  equation: 


C  V  -  C  u  = 
z  T  r  T 


ep 


-  C  XV  +  C  XU  +  — (C  Ti~Cn)P  3. 1.3.1 
z  n  r  n  p  rz  ^z 'r 


Equation  3. 1.3.1  involves  only  T-  and  p-  derivatives  of  the  state  variables 
and  may  be  marched  forward  according  to  the  same  prescription  as  we  described 
in  the  preceding  section  on  interior  points.  Then  the  application  of  the 
boundary  condition,  which  requires  the  normal  velocity  component  to  vanish, 
yields  updated  values  of  u  and  v,  A  similar  procedure  is  followed  for  the 
solid  phase  and,  for  both  phases,  on  an  p  =  constant  boundary. 

Then,  the  rate  of  propagation  of  granular  disturbances  is  computed.  If 
it  is  zero,  the  porosity  is  updated  directly  from  the  solid  phase  continuity 
equation  using  one  sided  differences  for  the  normal  derivatives  as  prescribed 
by  3. 1.2.4  and  3. 1.2.5.  If  the  rate  of  propagation  is  not  zero,  we  use  equa¬ 
tion  2.4.4.20  to  integrate  a  and  then  e  follows  from  2.4.4.13. 

In  any  case,  the  pressure  p  follows  from  equation  2.4.4.21  and  the  den¬ 
sity  may  then  be  determined  from  equation  2.4.4.12. 

In  discretizing  the  characteristic  forms  we  use  exactly  the  prescription 
of  the  previous  section,  on  interior  points,  to  evaluate  all  derivatives 
along  the  boundary.  Except  as  noted  in  section  3.1.5,  the  normal  spacewise 
derivatives  are  evaluated  according  to  3. 1.2.4  and  3. 1.2.5  with  the  sign 
chosen  to  ensure  differencing  inside  the  computational  domain.  The  coordin¬ 
ate  transformation  derivatives  are  deduced  using  centered  differences  for 
the  derivatives  along  the  boundary  and  the  usual  second  order  one  sided  three 
point  difference  for  the  normal  derivatives. 

3.1.4  Integration  at  Corner  Points 

At  the  corners  both  the  normal  derivatives  and  the  derivatives  along  the 
boundaries  must  be  resolved  in  one-sided  forms.  Moreover,  the  definition  of 
the  normal  direction  is  ambiguous  as  it  may  coincide  with  either  of  the  in¬ 
tersecting  coordinate  lines.  In  the  present  study  we  have  chosen  the  charac¬ 
teristic  plane  according  to  the  structure  of  the  flow  near  the  boundary.  The 
normal  direction  is  chosen  as  either  ^  -  r  or  as  p  -  x  according  as: 


k  u  +  i;  V 
'  z  5  r  5  ^ 


In  u  +  n  V  I 
'  z  p  r  p  ' 


Except  as  noted  in  the  following  section,  equations  3. 1.2.4  and  3. 1.2.5 
are  used  to  discretize  both  the  normal  and  the  tangential  spacewise  deriva¬ 
tives  of  the  state  variables.  The  mesh  transformation  derivatives  are  de¬ 
duced  according  to  the  usual  second  order  accurate  one  sided  difference 
formula . 
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3.1.5  Numerical  Devices 


The  purpose  of  this  section  is  to  record  certain  details  of  the  present 
method  of  solution  which  amount  to  ad  hoc  devices.  Two  topics  are  addressed. 
The  first  relates  to  the  ignition  of  particles  on  the  boundaries  and  the 
second  relates  to  the  representation  of  normal  derivatives  in  the  character¬ 
istic  analysis. 

1 8 

We  have  noted  previously  in  one  dimensional  studies  that  ignition 
can  be  significantly  delayed  for  particles  on  impermeable  boundaries  where 
the  slip  velocity  and  hence  the  convective  heating  may  vanish.  Physically, 
flamespreading  is  continued  to  the  boundary  by  different  mechanisms  than 
that  expressed  by  our  model.  In  our  one  dimensional  work  we  simply  extrapo¬ 
lated  the  flame  trajectory  to  the  boundary  in  order  to  estimate  the  ignition 
delay  for  particles  on  the  boundary.  In  the  present  two  dimensional  study, 
an  extrapolation  of  the  path  of  flamespreading  is  difficult  due  to  the  geo¬ 
metrical  factors.  Moreover,  in  subsequent  applications  to  bag  charges  we 
expect  to  deal  v/ith  permeable  boundaries  for  which  the  extended  ignition 
delay  does  not  arise.  Accordingly,  we  have  simply  compared  the  surface  tem¬ 
perature  of  a  particle  on  the  boundary  with  that  of  its  neighbor  in  the  in¬ 
terior  and  replaced  it  by  the  greater  of  the  two  values.  Thus  if  there  is 
significant  convective  heating  due  to  a  tangential  flow,  ignition  of  a 
boundary  particle  will  follow  in  the  natural  manner.  On  the  other  hand,  if 
the  boundary  corresponds  to  a  stagnation  region,  the  particle  temperature 
follows  that  of  its  neighbor  in  the  interior. 

As  described  in  section  3.1.3,  the  corrector  step  in  the  analysis  of  the 
boundary  values  involves  a  three  point  difference  formula  to  provide  overall 
second  order  formal  accuracy.  We  have  found  that  as  the  flamefront  approaches 
the  boundary,  the  three  point  formula  provides  less  real  accuracy  than  a  uni¬ 
formly  first  order  scheme.  The  same  effect  has  been  noted  in  our  one  dimen¬ 
sional  solutionslS.  Similar  observations  have  been  made  in  the  context  of 
highly  structured  single  phase  flow^*^.  Accordingly,  we  only  use  the  three 
point  corrector  when  ignition  has  occurred  at  the  boundary  point  in  question. 
During  the  pre-heat  period  we  use  the  first  order  difference  formula  on  both 
the  predictor  and  the  corrector  steps. 

3.2  Numerical  Results 


We  now  present  numerical  results  for  three  problems.  In  the  first  two 
problems  we  assume  that  the  projectile  does  not  intrude  into  the  combustion 
chamber.  These  serve  as  a  baseline  against  which  to  consider  the  third  prob¬ 
lem  in  which  the  projectile  is  assumed  to  have  a  long  boattail.  The  first 
two  problems  differ  from  one  another  only  in  respect  to  the  law  governing 
the  intergranular  stresses.  In  the  first  problem  the  stress  is  assumed  to 
be  reversible  while,  in  the  second,  it  is  taken  to  be  irreversible.  The 
data  bases  for  these  solutions  are  presented  in  section  3.2.1.  The  solu¬ 
tions  themselves  are  respectively  discussed  in  sections  3.2.2  through  3.2.4. 


■  kbboXt,  M.J.  "BoundoAy  CondUXdon  Calculation  Paoccdu/ici  Invi&cid 
SupeAtontc  Flow  Field-i"  Paoc.  hit  AIAA  Comp.  Fluid  Vynamict  Con^.  7973 
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3,2.1  Discussion  of  Data  Bases 


Figure  3.2. 1.1  presents  the  boundary  configurations  and  the  initial  mesh 
distributions  for  the  zero  intrusion  projectile,  discussed  in  sections  3.2.2 
and  3.2.3,  and  for  the  long  intrusion  projectile,  discussed  in  section  3.2.4. 
In  both  cases  the  mesh  has  been  generated  using  Dirichlet  data  along  each 
boundary  element. 

The  mesh  for  the  zero  intrusion  projectile  involves  16  points  in  the 
axial  direction  and  7  in  the  radial  direction.  An  overrelaxation  coefficient 
of  1.6  was  used  to  produce  a  mesh  convergent  to  1  part  in  10^  after  36  iter¬ 
ations.  The  21x7  mesh  used  for  the  long  intrusion  projectile  required  39 
iterations  to  converge  to  the  same  precision.  Notable  features  of  the 
physical  domain,  in  both  cases,  are  the  representation  of  the  curvature  of 
the  breech  and  the  taper  of  the  tube.  It  is  important  to  note  that  only  the 
four  corners  defined  by  the  intersections  of  the  four  boundary  elements  are 
treated  explicitly  in  the  calculation.  Thus,  in  the  case  of  the  long  intru¬ 
sion  projectile  the  corner  defined  by  the  termination  of  the  boattail  at  the 
base  is  perceived  implicitly  by  the  numerical  method  as  a  point  on  a  smooth 
and  continuously  differentiable  boundary  curve.  Therefore,  numerical  diffu¬ 
sion  is  to  be  expected  at  this  corner  with  a  corresponding  degradation  of 
accuracy . 

The  thermophysical  data  used  to  determine  the  three  solutions  are  pre¬ 
sented  in  Table  3.2, 1.1.  The  tabular  representation  of  the  igniter  discharge 
is  given  in  Table  3.2. 1.2.  Tables  3.2. 1,3  and  3.2. 1.4  respectively  present 
the  tabular  data  used  to  define  the  boundary  elements  for  the  zero  intrusion 
and  long  intrusion  cases. 

It  should  be  noted  that  the  data  bases  are  nominal  in  the  sense  that  we 
do  not  attempt  to  simulate  a  specific  in-service  round.  The  purpose  of  these 
calculations  is  simply  to  enable  us  to  appraise  the  operability  of  the  numer¬ 
ical  scheme . 

Attention  is  drawn  to  Table  3. 2. 1.2.  The  distribution  of  the  igniter 
venting  function  is  seen  to  be  quite  localized,  especially  in  the  radial  di¬ 
rection.  As  only  seven  points  are  used  to  resolve  the  radial  structure  of 
the  flow  we  should  not  anticipate  great  accuracy  in  these  solutions.  How¬ 
ever,  problems  of  resolution  concerning  the  igniter  venting  characteristics 
are  not  a  concern  in  the  present  study.  We  intend,  in  subsequent  work,  to 
pose  more  properly  the  physical  model  of  the  center  core  igniter.  Ultimately, 
it  will  be  treated  as  a  separate  computational  zone,  linked  to  the  behaviour 
of  the  two  phase  flow  in  the  propelling  charge  by  an  explicit  representation 
of  the  boundary  separating  the  two  zones. 
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Table  3. 2. 1.1  Thermo physical  Properties  Used  in  Two  Dimensional  Simulations 


Initial  Temperature  (°K)  305.6 

Initial  Pressure  (MPa)  0.1014 

Initial  Porosity  (-)  0.5 

Settling  Porosity  of  Granular  Bed  (-)  0.5 

Speed  of  Granular  Compression  Wave  (m/sec)  442. 

Speed  of  Granular  Expansion  Wave  (m/sec)*  1270. 

Density  of  Solid  Phase  (gm/cc)  1.6508 

Thermal  Conductivity  (J/cm-sec-^K)  0,0016 

Thermal  Diffusivity  (cm^/sec)  0.0006 

Ratio  of  Specific  Heats  (-)  1.25 

Molecular  weight  (gm/gmol)  28.8 

Covolume  (cc/gm)  1.084 

Ignition  Temperature  (°K)  488.9 

Chemical  Energy  of  Propellant  (J/gm)  4982. 

Burn  Rate  Additive  Constant  (cm/sec)  0. 

Burn  Rate  Pre-Exponential  Factor  (cm/sec-MPa  )  0.25159 

Burn  Rate  Exponent  (-)  0.6 

External  Diameter  of  Grain  (cm)  1.270 

Length  of  Grain  (cm)  2.54 

Diameter  of  Perforations  (cm)  0.127 

Number  of  Perforations  (-)  7 

Energy  of  Igniter  (J/gm)  3487 


Set  equal  to  442  m/sec  for  problem  with  reversible  stress  law 


/ 
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Table  3. 2. 1.2  Tabular  Data  Used  to  Define  Igniter  Discharge* 


Rate  of  Discharge  at  0  msec 


Radial  Location  (cm) 
Axial  Location  (cm) 

0. 

10. 

20. 


0 


1.5 


27.68 

27.68 

0 


27.68 

27.68 

0 


4.0 


0 

0 

0 


Rate  of  Discharge  at  1  msec 


Radial  Location  (cm) 
Axial  Location  (cm) 

0. 

10. 

20. 


0 


1.5 


27.68 

27.68 

0 


27.68 

27.68 

0 


4.0 


0 

0 

0 


Values  linearly  interpolated  for  arguments  within  table  range. 
Discharge  zero  outside  table  argument  range. 


Table  3.2. 1.3  Tabular  Data  Used  to  Define  Geometry  of  Chamber  with 
Zero  Intrusion  Projectile* 

Axial  Position  (cms)  Radial  Position  (cms) 


Breech 

3.0 

3.0 

2.5 

1.5 

0. 


0 

4.5 
6.0 

7.5 
9.0 


Projectile  Base 

33.0  0. 

33.0  8. 


Internal  Boundary 

3.0  0. 

33.0  0. 


External  Boundary 

0.  9.0 

33.0  8.0 


Points  located  on  each  boundary  element  by  linear  interpolation. 
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Table  3.2. 1.4  Tabular  Data  Used  to  Define  Geometry  of  Chamber  with 
Long  Intrusion  Projectile* 

Axial  Position  (cms)  Radial  Position  (cms) 

Screech 

0. 

4.5 

6,0 

7.5 
9.0 


Projectile  Base 

0. 

4.0 

8.0 


Internal  Boundary 

3.0  0. 

33,0  0. 


External  Boundary 

0.  9.0 

45.0  8.0 


Points  located  on  each  boundary  element  by  linear  interpolation. 


33.0 

33.0 

45.0 


3.0 

3.0 

2.5 

1.5 

0. 
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3.2.2  Zero  Intrusion  Projectile  with  Reversible  Granular  Stress  Law 


Some  aspects  of  this  solution  are  depicted  in  figures  3.2.2. 1  through 
3.2. 2.7.  These  are  confined  to  some  distributions  of  porosity  and  pressure 
at  three  times  during  the  flamespreading  process  and  the  distribution  of  gran¬ 
ular  stress  at  the  conclusion  of  flamespreading. 

From  the  distribution  of  igniter  venting  given  in  Table  3.2. 1.2  we  antici¬ 
pate  the  following  sequence  of  physical  events.  As  the  rate  of  venting  is 
quite  vigorous,  that  part  of  the  bed  adjacent  to  the  region  of  venting  is  sub¬ 
jected  to  intense  convective  heating  from  the  initial  instant.  Therefore, 
ignition  occurs  very  quickly  around  the  venting  region.  The  convective  flame 
sweeps  rapidly  outward  in  the  radial  direction  and  induces  ignition  of  the  en¬ 
tire  rear  portion  of  the  charge.  The  equilibration  of  pressure  over  the  cross- 
section  of  the  tube  occurs  faster  than  the  forward  spreading  of  the  flame. 

The  flow  soon  develops  a  one  dimensional  character,  with  the  distribution  of 
pressure  resembling  that  of  the  NOVA  code. 

Figures  3.2.2. 1,  3. 2.2, 2  and  3. 2. 2.3  present  the  distributions  of  poros¬ 
ity  at  three  times.  Flamespreading  is  complete  by  1.0  msec,  the  latest  time 
considered  in  these  figures. 

We  note  particularly  the  progressive  rarefaction  of  the  solid  phase  in 
the  region  of  igniter  blast  due  to  the  gas  dynamic  forces  exerted  by  the  expand¬ 
ing  products  of  combustion.  This  rarefaction  is,  of  course,  balanced  by  a  net 
compression  in  the  outer  regions.  The  degree  of  compaction  is  comparatively 
mild,  due  to  the  influence  of  the  granular  stress.  In  fact,  the  present  prob¬ 
lem  is  made  relatively  difficult  from  a  computational  point  of  view  by  the  ab¬ 
sence  of  an  explicit  treatment  of  ullage.  The  problem  would  be  better  posed, 
computationally  as  well  as  physically,  if  the  igniter  were  taken  to  vent  into 
the  bed  through  a  permeable  boundary.  In  such  a  case,  the  region  of  rarefac¬ 
tion  seen  in  figures  3.2.2. 1  through  3, 2.2.3  would  be  replaced  by  a  separate 
computational  zone  and  linked  to  the  compacted  bed  by  an  explicitly  repre¬ 
sented  internal  boundary.  This  aspect  of  the  solution  will  be  present  in  all 
the  three  cases  which  we  consider.  We  will,  in  the  succeeding  sections,  look 
a  little  more  closely  at  the  evidence  of  numerical  strain  imposed  by  the  re¬ 
quirement  of  capturing,  with  a  coarse  mesh,  the  interaction  between  the  com¬ 
pression  wave  and  the  rarefaction  wave  in  the  solid  phase. 

Concerning  the  figures,  we  make  the  following  general  observations.  They 
are  all  oblique  views  of  the  three  dimensional  surface  formed  by  the  state 
variable  in  question  and  the  axial  and  radial  coordinates.  The  axially  di¬ 
rected  lines  are,  in  fact,  contours  of  constant  r\  while  the  radially  directed 
lines  constitute  contours  of  constant  Here  ^  and  f]  are  the  computational 
coordinates.  The  intersections  of  the  two  families  define  the  mesh  point 
values  of  the  quantity  in  question.  The  contours  are  all  drawn  as  straight 
line  segments  between  successive  pairs  of  mesh  points  and  hidden  lines  are  re¬ 
moved  from  the  drawings.  As  the  computational  coordinates  follow  the  motion 
of  the  solid  phase,  the  mesh  distortion  can  be  inferred  from  these  figures, 
particularly  those  for  which  the  variable  in  question  is  essentially  uniform. 

It  should  also  be  noted  that  the  axial  length  scale  is  foreshortened  relative 
to  the  radial  length  scale  in  all  the  figures. 
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Returning  now  to  the  specific  results,  we  consider  figures  3. 2. 2. 4 
through  3.2.2. 7.  We  note  the  fully  two  dimensional  nature  of  the  flow  at 
0.2  msec  as  the  flame  is  still  in  the  process  of  spreading  radially.  By 
0.4  msec,  radial  equilibration  is  nearly  complete  and  by  1.0  msec,  the  pres¬ 
sure  distribution  has  virtually  no  radial  structure. 

Figure  3. 2. 2. 7  presents  the  distribution  of  granular  stress  at  the  con¬ 
clusion  of  flamespreading.  It,  too,  shows  very  little  radial  structure  at 
this  time. 

We  will  explore  more  fully  the  radial  structure  of  the  remaining  state 
variables  in  the  two  succeeding  sections.  Physically  speaking,  we  are  more 
interested  in  the  cases  in  which  the  granular  stress  law  is  irreversible.  We 
have  now  sufficient  information  about  the  present  problem  for  our  purpose, 
namely  to  provide  a  benchmark  comparison  with  the  next  problem. 

3.2.3  Zero  Intrusion  Projectile  with  Irreversible  Granular  Stress  Law 


The  solution  for  this  problem  is  represented  by  figures  3.2.3. 1  through 
3.2.3.33  and  is  therefore  displayed  in  ranch  greater  detail  than  that  of  the 
preceding  section. 

We  emphasize  that  the  only  physical  difference  between  the  present  prob¬ 
lem  and  that  of  the  preceding  section  relates  to  the  granular  stress  law.  In 
the  preceding  section  the  unloading  modulus  was  taken  to  be  the  same  as  the 
loading  modulus.  In  the  present  section,  as  shown  by  Table  3.2. 1.1,  the  un¬ 
loading  modulus,  proportional  to  the  square  of  the  wave  speed,  is  approxi¬ 
mately  ten  times  as  great  as  the  loading  modulus. 

Figure  3.2.3. 1  presents  the  distribution  of  porosity  following  the  com¬ 
pletion  of  flamespreading  at  1.0  msec.  By  comparison  with  figure  3.2. 2.3  we 
see  that  the  relative  absence  of  elastic  recovery  in  the  present  example  re¬ 
sults  in  somewhat  greater  rarefaction  in  the  region  directly  influenced  by 
the  primer  blast.  Figure  3.2.3. 2  presents  the  porosity  at  a  still  later 
time  and  shows  even  greater  expansion.  Indeed,  the  porosity  is  very  nearly 
discontinuous  at  the  leading  edge  of  the  granular  rarefaction. 

Figures  3. 2.3.3  through  3. 2. 3. 6  present  distributions  of  pressure  at 
various  times.  The  first  three  of  these  may  be  compared  with  their  counter¬ 
parts  in  the  previous  section.  Evidently,  the  change  in  the  granular  stress 
law  has  produced  no  qualitative  changes  in  the  pressure  distributions  and 
only  minor  quantitative  differences  are  seen.  Figure  3. 2.3.6  shows  that  by 
1.3  msec  the  pressure  has  completely  equilibrated  throughout  the  chamber  in 
spite  of  the  transient  flamespreading  process.  This  may  be  regarded  as  an 
indication  of  the  high  degree  of  mechanical  dissipation  present  in  the  two 
phase  flow.  The  figure  also  enables  one  to  assess  the  total  mesh  distortion 
which  takes  place  during  flamespreading. 

Thus  far  we  have  discussed  flamespreading  without  providing  any  quanti¬ 
tative  details  of  its  history.  In  figures  3. 2.3.7  through  3,2.3.11  we  pre¬ 
sent  the  distributions  of  the  surface  temperature  of  the  solid  phase.  The 
fully  two  dimensional  flamespreading  induced  by  the  igniter  blast  is  seen 
clearly  in  figures  3.2.3. 7  and  3.2. 3.8.  By  0.3  msec,  as  shown  in  figure 
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3. 2. 3, 9,  flamespreading  has  a  largely  one  dimensional  character  although  the 
continuing  venting  of  the  igniter  produces  a  somewhat  greater  rate  of  spread¬ 
ing  at  the  centerline.  This  is  also  seen  in  figure  3.2.3.10.  By  1.0  msec, 
as  seen  in  figure  3.2.3.11,  flamespreading  is  complete. 

These  figures  indicate  just  how  short  the  preheat  or  induction  phase  is 
at  the  leading  edge  of  the  flame.  The  surface  temperature  distribution  es¬ 
sentially  consists  of  two  plateaux,  the  lower  at  ambient  temperature  and  the 
higher  at  the  ignition  temperature.  Actually,  the  upper  surface  is  seen  to 
depart  slightly  from  an  ideally  uniform  distribution.  The  departure  is  due 
to  a  slight  overshoot  of  the  temperature  history  during  the  time  step  corres¬ 
ponding  to  the  ignition  event  at  any  given  mesh  point. 

In  figures  3.2.3.12  through  3.2.3.15,  we  consider  the  distributions  of 
gas  density  at  four  times.  These  exhibit  a  pronounced  entropy  layer  near  the 
external  and  forward  boundaries  when  the  pressure  distributions  are  also  taken 
into  account.  The  density  at  the  external  and  forward  boundaries,  where  igni¬ 
tion  is  delayed  due  to  weak  convection,  is  seen  to  be  very  much  greater  than 
that  elsewhere  in  the  chamber.  In  effect,  the  density  at  these  boundaries  re¬ 
sults  largely  from  an  isentropic  response  to  the  pressure  field.  The  tempera¬ 
ture  distribution  is  correspondingly  depressed,  as  shown  in  figure  3.2.3.16. 

Of  course,  once  ignition  has  occurred  on  the  boundaries  the  density  distribu¬ 
tion  tends  to  become  more  uniform  as  may  be  seen  by  comparing  figures  3.2.3.14 
and  3.2.3.15. 

Figures  3.2.3.12  and  3.2.3.13  reveal  the  rather  weak  adiabatic  precursor 
that  runs  ahead  of  the  convective  flame  in  the  axial  direction.  It  is  also 
clear  why  upstream  differencing  of  the  convective  terms  is  necessary  for 
stability.  The  use  of  the  standard  MacCormack  prescription  for  the  convective 
terms  is  found  to  lead  to  negative  values  of  the  density  at  mesh  points  adja¬ 
cent  to  the  cool  boundaries. 

We  turn  now  to  some  additional  aspects  of  the  behaviour  of  the  solid 
phase,  looking  first  at  some  distributions  of  granular  stress  and  subsequently 
at  the  velocity  components. 

The  distributions  of  granular  stress  are  presented  in  figures  3.2.3.17 
through  3.2.3.22.  In  all  of  these  figures  the  granular  stress  is  seen  to  van¬ 
ish  in  the  region  of  igniter  venting  where  the  porosity  corresponds  to  dispersed 
flow.  Figures  3.2.3.17  through  3.2.3.19  indicate  that  quite  a  strong  radial 
compression  is  established  by  the  igniter  blast.  Then,  as  the  flame  penetrates 
to  the  outside  of  the  charge,  the  granular  stress  quickly  unloads.  By  0.6  msec, 
as  shoxNm  in  figure  3.2.3.20,  we  have  an  essentially  one  dimensional  compression 
wave  moving  in  the  axial  direction.  The  granular  stress  exerted  on  the  base 
of  the  projectile  is  seen,  in  figure  3.2.3.21,  eventually  to  rise  to  a  value 
comparable  to  that  experienced,  much  earlier,  at  the  external  boundary  near 
the  breech  in  figure  3.2.3.19.  By  1.3  msec,  figure  3.2.3.22,  the  granular 
stress  has  decayed  strongly  in  the  forward  region.  At  this  point  the  grains 
have  rebounded  somewhat  and  a  mild  level  of  stress  is  seen  in  the  outer  part 
of  the  breech. 

In  figure  3.2.3.22  we  see,  for  the  first  time,  an  unsatisfactory  level  of 
numerical  noise.  This  is  due  in  part  to  the  non-analytical  character  of  the 
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granular  stress  law  and  also  to  the  numerical  strain  imposed  by  the  interac¬ 
tion  of  the  reflected  compression  wave  with  the  rarefaction  in  the  region  of 
igniter  blast.  We  do  not  address  the  problem  of  alleviating  these  nijmerical 
wiggles  in  the  present  study.  In  the  first  place,  they  are  associated  with 
relatively  low  amplitudes  of  granular  stress.  In  the  second  place,  as  we 
have  already  noted,  subsequent  studies  will  introduce  an  explicit  internal 
discontinuity  to  separate  the  region  of  igniter  venting  from  the  compacted  bed. 

Some  idea  of  the  granular  velocity  field  may  be  had  by  reference  to 
figures  3.2.3.23  through  3.2.3.26.  The  maximum  grain  velocity  appears  to  oc¬ 
cur  at  the  periphery  of  the  region  of  igniter  venting.  We  see  that  both  the 
axial  and  radial  components  have  maxima  equal  to  11  m/sec  and  that  these 
maxima  are  virtually  constant  throughout  flamespreading.  Virtually  no  radial 
motion  of  the  charge  occurs  forward  of  the  region  of  igniter  venting.  The 
axial  component  is  seen  not  to  exceed  5  m/sec  in  this  region  at  0.6  msec  and, 
by  1.0  msec,  to  be  slightly  negative  as  the  axial  compression  wave  is  reflec¬ 
ted  from  the  base  of  the  projectile. 

We  conclude  our  discussion  of  this  solution  by  examining  some  distribu¬ 
tions  of  the  components  of  the  gas  velocity.  The  axial  component  is  illus¬ 
trated  at  two  times  in  figures  3.2.3.27  and  3.2.3.28.  These  reveal  no  espec¬ 
ially  surprising  features  in  view  of  the  preceding  discussion.  The  gas  vel¬ 
ocity  is  seen  to  have  a  value  of  approximately  220  m/sec  at  0.6  m.sec.  How¬ 
ever,  the  maximum  value  has  subsided  to  approximately  110  m/sec  by  1.0  msec. 

The  largest  value  of  the  axial  component  occurs  shortly  after  the  start  of 
igniter  venting  and  is  approximately  270  m/sec  or  slightly  more  than  that 
shown  in  figure  3.2.3.27. 

The  distributions  of  the  radial  component  of  the  gas  phase  velocity  are 
of  somewhat  greater  interest,  figures  3.2.3.29  through  3.2.3,33.  At  0,2  msec, 
the  maximum  value  of  the  radial  component  is  approximately  270  m/sec  which  is 
the  same  as  that  for  the  axial  component.  By  0.6  msec,  the  maximum  has  al¬ 
ready  decreased  to  about  80  m/sec,  by  1.0  msec  to  about  30  m/sec  and  by  1.3 
msec  to  less  than  10  m/sec.  The  rapid  damping  of  the  radial  component  of  the 
gas  velocity  is  in  accord  with  our  previous  findings  based  on  the  solutions 
for  cylindrical  flow^  . 

It  is  interesting  to  note,  in  figures  3.2.3.30  and  3.2.3.31,  the  exis¬ 
tence  of  two  maxima  in  the  radial  velocity  distributions.  The  profile  near 
the  breech  is  associated  with  the  venting  of  the  igniter  which  is  represented 
as  continuing  for  1  msec.  The  forward  profile  is  apparently  driven  by  the 
tendency  of  the  convective  flame  to  produce  ignition  somewhat  earlier  at  the 
centerline  than  at  the  outside.  Once  the  igniter  venting  is  terminated  the 
radial  velocity  distribution  at  the  rear  reverses  itself.  Gas  now  flows  in¬ 
ward  to  supply  the  region  attenuated  by  the  igniter  blast  and  which  is  accord¬ 
ingly  deficient  in  terms  of  local  gas  generation.  This  is  seen  in  figures 
3.2.3.32  and  3.2.3.33,  The  forward  profile  disappears  altogether  once  igni¬ 
tion  of  the  bed  is  complete. 
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3.2.4  Long  Intrusion  Projectile 


This  problem  differs  from  that  of  the  preceding  section  in  respect  to 
the  geometry  of  the  projectile  which  is  now  taken  to  have  a  boattail  which 
intrudes  into  the  combustion  chamber.  Thus  we  introduce  an  added  measure  of 
geometrical  complexity . 

In  many  respects,  however,  the  solution  for  the  present  problem  resembles 
that  of  the  preceding  section.  Therefore,  we  will  not  present  as  many  details 
of  the  solution  as  we  did  in  the  preceding  case.  We  will  first  provide  some 
indications  of  the  similarity  of  the  two  solutions  and,  subsequently,  dwell 
on  the  differences.  The  present  solution  is  represented  by  figures  3. 2.4.1 
through  3.2.4.21. 

Figure  3.2.4. 1  presents  the  distribution  of  porosity  at  1.4  msec,  follow¬ 
ing  the  conclusion  of  flamespreading.  The  additional  delay  to  complete  flame¬ 
spreading  is,  of  course,  due  to  the  extension  of  the  chamber  defined  by  the 
region  around  the  boattail.  We  have,  as  in  the  preceding  cases,  a  pronounced 
rarefaction  in  the  region  of  igniter  blast  accompanied  by  a  modest  compression 
of  the  remainder  of  the  charge. 

Figures  3. 2.4. 2  through  3. 2. 4. 5  present  the  distributions  of  pressure  at 
four  different  times.  These  resemble,  qualitatively,  the  results  presented 
in  the  previous  section.  In  particular  we  see  that  shortly  after  the  conclu¬ 
sion  of  flamespreading,  figure  3. 2.4.5,  the  pressure  is  virtually  uniform 
throughout  the  chamber,  including  the  region  around  the  boattail.  Figure 
3. 2.4.5  also  provides  a  clear  indication  of  the  total  mesh  distortion  which 
occurs  during  flamespreading. 

We  note  that  there  is  no  sign  of  numerical  wiggles  in  the  vicinity  of  the 
corner  defined  by  the  intersection  of  the  boattail  and  the  base  of  the  projec¬ 
tile,  The  density  and  temperature  of  the  gas  are  equally  well  behaved  as  in¬ 
dicated  by  figures  3. 2. 4. 6  and  3. 2. 4. 7.  Of  course,  the  pronounced  entropy 
layer  is  still  present  although,  due  to  the  somewhat  later  time  considered 
here,  the  effect  is  diminished  at  the  external  boundary. 

Figures  3. 2.4. 8  through  3.2.4.11  illustrate  the  path  of  flamespreading. 

By  1.0  msec  the  entire  portion  of  the  charge  to  the  rear  of  the  base  is  fully 
ignited,  in  keeping  with  the  finding  of  the  previous  section. 

Distributions  of  granular  stress  are  shown,  at  three  times,  in  figures 
3.2.4.12  through  3.2.4.14.  The  distribution  at  0.3  msec  is  seen  to  be  vir¬ 
tually  the  same  as  in  the  previous  case,  as  one  would  expect.  By  1.0  msec 
some  differences  are  apparent  as  the  stress  becomes  quite  large  in  the  region 
around  the  boattail.  At  1.3  msec  the  stress  in  the  corner  defined  by  the  in¬ 
tersection  of  the  tube  and  the  boattail  reaches  its  maximum  value  of  nearly 
20  MPa.  Subsequently,  unloading  occurs,  partly  as  a  consequence  of  wave  re¬ 
flection  and  partly  due  to  combustion  of  the  grains.  We  note  that  20  MPa 
represents  a  rather  large  value  of  the  average  granular  stress.  Considering 
that  the  actual  bearing  surface  to  support  the  stress  may  be  quite  small,  the 
localized  values  within  the  grains  may  be  significantly  higher  than  this  value. 
Therefore  the  possibility  of  grain  fracture  ought  not  be  overlooked. 
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We  conclude  by  examining  some  of  the  distributions  of  the  velocity  com¬ 
ponents  of  each  of  the  phases.  The  axial  component  of  the  gas  velocity  is 
shown  in  figures  3.2.4.15  through  3.2.4.17.  Again,  the  distribution  at  0.6 
msec  resembles  closely  that  of  the  previous  section,  the  gas  being  quiescent 
in  the  region  around  the  boattail  at  this  time.  By  0.1  msec,  the  convection 
around  the  boattail  has  strengthened  as  the  flame  begins  to  penetrate  that 
region.  By  1.4  msec,  however,  the  axial  component  of  the  gas  phase  velocity 
has  subsided  significantly  and  a  weak  backflow  is  apparently  about  to  occur, 
near  the  breech,  in  the  region  rarefied  by  the  igniter  blast, 

The  corresponding  distributions  of  the  radial  component  of  the  gas  phase 
velocity  are  shown  in  figures  3.2.4.18  through  3.2.4.20.  The  only  noteworthy 
comment  concerns  the  existence  of  a  numerical  wiggle  in  the  corner  defined  by 
the  tube  and  the  boattail.  This  is  probably  due  to  the  relatively  coarse  mesh 
spacing  in  that  region,  compounded  by  the  delayed  ignition.  Once  flamespread¬ 
ing  is  complete,  the  wiggle  is  seen  to  be  much  less  pronounced,  figure  3.2.4.20. 

Finally,  in  figures  3.2.4.21  and  3.2.4.22  we  present  the  distributions  of 
the  axial  and  radial  components  of  the  solid  phase  velocity.  These  have  the 
expected  correspondence  with  their  counterparts  of  the  preceding  section. 
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4.0  CONCLUDING  REMARKS 


We  have  described  the  overall  theoretical  framework  for  a  model  of  the 
two  dimensional  flow  in  a  gun.  The  objective  of  the  model  being  the  analysis 
of  longitudinal  wave  propagation,  emphasis  has  been  given  to  the  precise  nu^ 
merical  treatment  of  such  potentially  important  aspects  of  the  propelling 
charge  as  the  location  of  ullage  and  of  the  influence  of  bag  material. 

As  a  step  towards  the  complete  numerical  implementation  of  the  model, 
we  have  encoded  a  method  of  solution  of  the  equations  of  two  dimensional  hetero 
geneous  reacting  two  phase  flow  and  presented  results  for  chamber  geometries 

whose  complexity  is  typical  of  what  we  expect  to  encounter  in  howitzer  charges. 

We  now  provide  some  comments  on  the  extendability  of  the  present  methods 
to  account  for  ullage  and  we  indicate  the  next  steps  in  the  path  of  develop¬ 
ment  of  the  code. 

Regarding  the  extension  of  the  present  method  to  account  for  ullage, 
we  comment  successively  on  the  overall  approach  to  programming  strategy,  the 
influence  of  the  gas  permeable  internal  boundaries  and  the  mapping  techniques 
for  the  regions  of  ullage. 

The  strategy  is  a  natural  extension  of  that  used  previously  in  the  devel¬ 
opment  of  the  one  dimensional  model^^.  We  assume  the  existence  of  two  user 
supplied  spatial  resolution  factors,  one  axial  and  one  radial.  Then  a  con¬ 
tinuum  representation  is  made  of  either  the  radial  or  the  axial  structure  of 
a  given  region  of  the  flow  in  accordance  with  a  criterion  which  compares  the 
extent  of  the  region  with  that  of  the  chamber  as  a  whole.  If  the  radial  ex¬ 
tent  of  the  region  exceeds  the  product  of  the  radial  resolution  factor  and 
the  radius  of  the  tube,  a  continuum  representation  is  made  of  the  radial 
structure  of  the  flow  and  likewise  for  the  axial  structure.  On  this  basis  a 
given  region  can  be  determined  as  two  dimensional,  quasi-one-dimensional ,  or  as 
lumped  parameter. 

The  present  method  is  well  suited  to  the  problem  of  the  gas  permeable 
internal  boundary.  The  jump  conditions  may  be  differentiated  along  the  path¬ 
line  of  the  boundary  to  provide  a  system  of  partial  differential  equations 
which  are  updated  explicitly  in  combination  with  the  characteristic  forms 
to  yield  new  values  of  the  state  variables  on  both  sides  of  the  internal 
boundary. 

As  regards  the  mapping  of  the  regions  of  ullage,  new  considerations 
arise  only  when  the  ullage  is  treated  as  a  two  dimensional  region.  However, 
the  equipotential  map  is  easily  made  time  dependent  so  that  no  difficulty  is 
anticipated  in  maintaining  a  suitable  mesh  distribution. 

We  conclude  with  some  comments  on  the  application  of  the  code,  in  its 
present  form,  to  the  analysis  of  flamespreading  in  a  charge  for  which  the 
projectile  intrusion  is  extremely  large.  Figure  4.1  presents  two  possible 
mesh  distributions  for  the  representation  of  a  105mm  Howitzer  with  the 
XM622  projectile.  The  first  distribution  was  generated  by  mapping  the  entire 
rear  surface  of  the  projectile,  including  afterbody  and  boattail  onto  one 
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side  of  the  square  computational  domain.  This  approach  is  seen  to  provide 
very  crowded  mesh  lines  in  the  rear  of  the  chamber  while  the  resolution  at 
the  front  is  poor.  Since  the  physical  problem  of  interest  relates  to  the 
possibility  of  large  granular  stress  in  the  corner  defined  by  the  intersec¬ 
tion  of  the  boattail  and  the  tube>  this  mesh  is  unsatisfactory. 

The  second  mesh  was  generated  by  mapping  only  the  boattail  onto  one 
side  of  the  computational  domain.  The  portion  of  the  centerline  between  the 
breech  and  the  base  of  the  projectile  together  with  the  entire  afterbody  has 
been  mapped  onto  the  computational  boundary  normally  thought  of  as  the  in¬ 
ternal  boundary.  This  is  seen  to  provide  a  much  better  distribution  of  mesh 
points.  As  motion  of  the  projectile  is  not  considered,  no  difficulty  arises 
from  the  apparent  splitting  of  the  physical  boundary  element. 

However,  it  should  be  noted  that  a  stagnation  boundary  condition  is 
applied  to  both  phases  at  the  corners  of  the  computational  mesh  and  nowhere 
else.  Accordingly,  further  attention  will  be  required  to  the  problem  of 
posing  correctly,  within  the  macroscopic  point  of  view,  the  boundary  condi¬ 
tions  at  sharp  internal  edges  of  the  type  associated  with  projectiles  having 
a  long  afterbody. 
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Fig.  3.2.2, 1  Distribution  of  porosity  at  0.2  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3.2.2, 2  Distribution  of  porosity  at  0.4  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3. 2. 2. 3  Distribution  of  porosity  at  1.0  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3. 2. 2.4  Distribution  of  pressure  at  0.2  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3. 2. 2. 5  Distribution  of  pressure  at  0.4  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3. 2. 2. 6  Distribution  of  pressure  at  1.0  msec  in  case  of  zero 
intrusion  projectile  with  reversible  granular  stress 
law 
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Fig.  3. 2. 3.1  Distribution  of  porosity  at  1.0  msec  in  case  of  zero 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3. 2. 3. 2  Distribution  of  porosity  at  1.3  msec  in  case  of  zero 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3. 2. 3. 3  Distribution  of  pressure  at  0.2  msec  in  case  of  zero 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3.2. 3.4  Distribution  of  pressure  at  0.4  msec  in  case  of  zero 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3. 2. 3. 5  Distribution  of  pressure  at  1.0  msec  in  case  of  zero 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3. 2. 3. 6  Distribution  of  pressure  at  1.3  msec  in  case  of 
intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3. 2. 3. 7  Distribution  of  solid  phase  surface  temperature  at 
0.1  msec  in  case  of  zero  intrusion  projectile  with 
irreversible  granular  stress  law 
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Fig.  3. 2. 3. 8  Distribution  of  solid  phase  surface  temperature  at 
0.2  msec  in  case  of  zero  intrusion  projectile  with 
irreversible  granular  stress  law 
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Fig.  3. 2. 3. 9  Distribution  of  solid  phase  surface  temperature  at 
0.3  msec  in  case  of  zero  intrusion  projectile  with 
irreversible  granular  stress  law 
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Fig.  3.2.3.10  Distribution  of  solid  phase  surface  temperature  at 
0.4  msec  in  case  of  zero  intrusion  projectile  with 
irreversible  granular  stress  law 
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Fig.  3.2.3.11  Distribution  of  solid  phase  surface  temperature  at 
1.0  msec  in  case  of  zero  intrusion  projectile  with 
irreversible  granular  stress  law 
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Fig.  3.2.3.12  Distribution  of  density  of  gas  at  0.2  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 


-  92 


DENSITT  (GM/CC)(X10-»  J 

0.01  O.OB  O.IS  0.11  0.20 


ZERO  INTRUSION  PROJECTILE 
IRREVERSIBLE  STRESS  LRN 
STEP  57 


TIHEmSED  O.UOO 


RXIRL  LOCRTION  (CHS) 


Fig.  3.2.3.13  Distribution  of  density  of  gas  at  0.4  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.14  Distribution  of  density  of  gas  at  1.0  msec  in  case 

of  zero  intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3.2.3.15  Distribution  of  density  of  gas  at  1.3  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.16  Distribution  of  temperature  of  gas  at  1.3  msec  in  case 
of  zero  intrusion  projectile  with  irreversible  granular 
stress  law 
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Fig.  3.2.3.17  Distribution  of  granular  stress  at  0.1  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.19  Distribution  of  granular  stress  at  0.3  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.20  Distribution  of  granular  stress  at  0.6  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.21  Distribution  of  granular  stress  at  1.0  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.22  Distribution  of  granular  stress  at  1.3  msec  in  case 
of  zero  intrusion  projectile  with  irreversible 
granular  stress  law 
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Fig.  3.2.3.23  Distribution  of  axial  component  of  solid  phase 
velocity  at  0.6  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 


VELOCITY  fCH/SECHX10> 


ZEAS  IHTRUSIBN  PABJECTILE 
IHREVEHSIBLE  STRESS  LAM 

STEP  IW  TIHEWSEC)  1.000 


Fig.  3.2.3.24  Distribution  of  axial  component  of  solid  phase 
velocity  at  1.0  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.25  Distribution  of  radial  component  of  solid  phase 
velocity  at  0.6  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.26  Distribution  of  radial  component  of  solid  phase 
velocity  at  1.0  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.27  Distribution  of  axial  component  of  gas  phase 

velocity  at  0.6  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.28  Distribution  of  axial  component  of  gas  phase 

velocity  at  1.0  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.29  Distribution  of  radial  component  of  gas  phase 
velocity  at  0.2  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.30  Distribution  of  radial  component  of  gas  phase 
velocity  at  0.6  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.31  Distribution  of  radial  component  of  gas  phase 
velocity  at  1.0  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.32  Distribution  of  radial  component  of  gas  phase 
velocity  at  1.1  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.3.33  Distribution  of  radial  component  of  gas  phase 
velocity  at  1.3  msec  in  case  of  zero  intrusion 
projectile  with  irreversible  granular  stress  law 
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Fig.  3.2.4. 1  Distribution  of  porosity  at  1.4  msec  in  case  of 
long  intrusion  projectile 
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Fig.  3. 2.4. 2  Distribution  of  pressure  at  0.2  msec  in  case 
of  long  intrusion  projectile 
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Fig.  3. 2.4. 3  Distribution  of  pressure  at  0.4  msec  in  case 
of  long  intrusion  projectile 
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Fig.  3. 2. 4. 4  Distribution  of  pressure  at  1.0  msec  in  case 
of  long  intrusion  projectile 
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Fig.  3. 2.4.5  Distribution  of  pressure  at  1.4  msec  in  case 
of  long  intrusion  projectile 
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Fig.  3. 2. 4. 6  Distribution  of  density  of  gas  at  1.4  msec  in  case 
of  long  intrusion  projectile 
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Fig.  3. 2. 4. 8  Distribution  of  solid  phase  surface  temperature 
at  0.2  msec  in  case  of  long  intrusion  projectile 
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Fig.  3. 2. 4. 9  Distribution  of  solid  phase  surface  temperature 
at  0.4  msec  in  case  of  long  intrusion  projectile 
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Fig.  3.2.4.10  Distribution  of  solid  phase  surface  temperature 
at  1.0  msec  in  case  of  long  intrusion  projectile 
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Fig.  3.2.4.11  Distribution  of  solid  phase  surface  temperature 
at  1.4  msec  in  case  of  long  intrusion  projectile 
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Fig.  3*2.4.12  Distribution  of  granular  stress  at  0.3  msec  in 
case  of  long  intrusion  projectile 
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Fig.  3.2.4.15  Distribution  of  axial  component  of  gas  phase 

velocity  at  0.6  msec  in  case  of  long  intrusion 
proj  ectile 
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Fig.  3.2.4.16  Distribution  of  axial  component  of  gas  phase 

velocity  at  1.0  msec  in  case  of  long  intrusion 
projectile 
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Fig.  3.2.4.17  Distribution  of  axial  component  of  gas  phase  velocity 
at  1.4  msec  in  case  of  long  intrusion  projectile 
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Fig.  3.2.4.18  Distribution  of  radial  component  of  gas  phase 
velocity  at  0.6  msec  in  case  of  long  intrusion 
proj  ectile 
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Fig.  3.2.4.19  Distribution  of  radial  component  of  gas  phase 
velocity  at  1.0  msec  in  case  of  long  intrusion 
projectile 
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Fig.  3.2.4.20  Distribution  of  radial  component  of  gas  phase 
velocity  at  1.4  msec  in  case  of  long  intrusion 
proj  ectile 
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Fig.  3.2.4.21  Distribution  of  axial  component  of  solid  phase 
velocity  at  1.0  msec  in  case  of  long  intrusion 
projectile 
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NOMENCLATURE 


Certain  symbols  are  used  in  different  contexts  with  different  meanings. 

It  is  not  expected  that  confusion  will  result  as  the  contexts  are  very  differ¬ 
ent  . 


mi,mo 


English  Symbols 


Cross  sectional  area  of  a  quasi-one-dimensional  flow 
Rate  of  propagation  of  granular  disturbances 
Value  of  a  for  settled  bed  during  compression 

Value  of  a  for  unloading  or  reloading  bed  when  porosity  is  less  than 

settling  porosity 

Burn  rate  additive  constant 

Burn  rate  pre-exponential  factor 

Covolume  of  gas  phase 

Speed  of  sound  in  gas  phase 

Specific  heats  at  constant  volume  and  constant  pressure 
Effective  diameter  of  a  grain  of  propellant 

Rate  of  surface  regression  of  burning  propellant 
Sum  of  internal  and  kinetic  energies 
Internal  energy  of  gas  phase 

Chemical  energy  released  in  combustion  of  solid  phase 
Bore  resistance 

Interphase  drag 

Constant  used  to  reconcile  units  of  measurement 

Parameter  used  to  deduce  propellant  surface  temperature  by  cubic 

profile  method 

Heat  transfer  coefficient 

Polar  moment  of  inertia  of  projectile 

Thermal  conductivity 

Projectile  Mass 

Effective  projectile  Mass 

Mass  production  per  unit  volume  per  unit  time  due  to  propellant 
combustion 

Mass  fluxes  to  and  from  a  region 

Normal  vector 
Burn  rate  exponent 
Prandtl  number 
Pressure 
Heat  flux 

Intergranular  stress 

Reynolds  number  based  on  effective  particle  diameter 

Radii  of  surfaces  of  quasi-one-dimensional  axial  flow  across  which 

mass  enters  and  exits,  respectively 

Radial  coordinate 

Surface  area  of  a  propellant  grain 
Surface  area  of  propellant  per  unit  volume 
Gas  temperature 

Surface  temperature  of  solid  phase 
Time 

Gas  velocity  vector,  components  (u,v) 
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Up  Solid  phase  velocity  vector,  components  (Up,Vp) 

Vp  Volume  of  a  propellant  grain 

w, Wp  C  "  component  of  gas,  solid  velocity  in  computational  plane 

x, Xp  ri  -  component  of  gas,  solid  velocity  in  computational  plane 

z  Axial  coordinate 

Greek  Symbols 

a  Characteristic  coordinate 

ap  Thermal  diffusivity  of  a  grain 

3  Virtual  mass  coefficient;  also  characteristic  coordinate 

Y  Ratio  of  specific  heats 

e  Porosity 

Eq  Settling  porosity 

^  Computational  coordinate 

r|  Computational  coordinate 

G  Angle  of  rifling  of  tube 

A  Coefficient  used  to  render  balance  equations  pseudo-totally 

hyperbolic 
y  Viscosity 

p  Density  of  gas 

Pp  Density  of  solid  propellant,  a  constant 

a  (l-e)R 

T  Time  coordinate  in  computational  frame 

ijj  Rate  of  production  of  gas  per  unit  volume  due  to  igniter;  also 

used  to  represent  a  column  vector  of  state  variables 

Special  S3mibols  and  subscripts 

D/Dt  Convective  derivative  along  average  gas  phase  streamline 

D/Dtp  Convective  derivative  along  average  solid  phase  streamline 

•  A  dot  over  a  quantity  indicates  differentiation  with  respect  to 

time  along  a  material  path  line 

IG,  p  The  subscript  IG  is  used  to  denote  properties  of  the  igniter  and 
p  is  used  to  denote  properties  of  the  solid  phase.  Gas  phase 
properties  are  unsubscripted 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet  and  return  it  to  Director,  US  Army  Ballistic  Research 
Laboratory,  ARRADCOM,  ATTN:  DRDAR-TSB,  Aberdeen  Proving  Ground, 
Maryland  21005.  Your  comments  will  provide  us  with  information 
for  improving  future  reports. 

1 .  BRL  Report  Number^ _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas ,  etc . )  _  _ _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5,  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 


Name: 

Telephone  Number: 
Organization  Address: 


